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Abstract 

We  consider  electromagnetic  interrogation  problems  for  complex  materials  involving  distributions  of  polarization 
mechanisms  and  also  distributions  for  the  parameters  in  these  mechanisms.  A  theoretical  and  computational  frame¬ 
work  for  such  problems  is  given.  Computational  results  for  specific  problems  with  multiple  Debye  mechanisms  are 
given  in  the  case  of  discrete,  uniform,  log-normal,  and  log-Bi-Gaussian  distributions. 

Keywords:  Electromagnetic  interrogation  with  pulsed  antenna  source  microwaves  and  inverse  problems,  complex 
dielectric  materials,  distributions  of  relaxation  parameters  and  mechanisms. 


1  Introduction 

For  at  least  the  past  century  [50,  51],  scientific  investigators  have  sought  to  understand  what  happens 
to  electromagnetic  fields  (and  how  to  mathematically  model  the  associated  phenomena)  when  they  are 
introduced  into  complex  materials  such  as  biotissue  and  more  general  dielectrics,  conductors  and  magnetics. 
More  specifically,  a  fundamental  question  is  how  to  model  dispersion  and  dissipation  of  the  fields  in  these 
complex  materials.  This  has  most  often  led  to  the  use  of  Maxwell’s  equations  in  a  non-vacuum  environment 
which  entails  constitutive  relationships  for  polarization  (in  dielectrics),  magnetization  (in  magnetic  materials) 
and  conductivity.  We  focus  here  on  modeling  polarization  in  dielectric  materials  for  which  we  develop  a  new 
modeling  framework.  Even  though  we  treat  only  polarization  as  our  dispersive  mechanism  in  our  formulation 
(adopting  Ohm’s  law  for  conductivity  and  considering  non-magnetic  materials),  the  approach  is  sufficiently 
general  so  as  to  be  readily  extended  to  treat  magnetization  and  conductivity  in  materials  (each  in  some  type 
of  convolution  representation  involving  susceptibility  kernels,  e.g.,  see  [2,  3]).  We  develop  a  framework  that 
allows  not  only  uncertainty  (through  distributions  of  parameters  representing  molecular  variability)  at  the 
molecular  level,  but  also  allows  for  the  presence  of  multiple  polarization  mechanisms  in  the  material. 

We  first  explain  a  conceptual  framework  in  the  context  of  the  3-D  Maxwell  equations  in  a  dielectric 
material.  In  particular,  we  use  Maxwell’s  equations  which  govern  the  electric  field  E  and  the  magnetic  field 
H  in  a  domain  T>  =  Qq  U  with  charge  density  p  in  the  material  fl  while  the  ambient  fio  is  treated  as  a 
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vacuum.  Thus  we  first  consider  the  system 


(i)  +  J  —  V  x  H  =  0,  in  (0,  T )  x  V, 

(ii)  —  +  VxE  =  0,  in  (0,T)  xD, 

(iii)  V  ■  D  =  p,  in  (0,  T)  x  2? , 

(iv)  V  •  B  =  0,  in  (0,T)  x  2?, 

(v)  E  x  n  =  0,  on  (0,  T)  x  <92?, 

,  (vi)  E(0,x)  =  0,  H(0,x)=0,  in  2?. 


As  usual,  the  current  J  is  composed  of  the  source  current  Js  and  the  conductive  current  Jc.  Within  the 
domain  we  have  constitutive  relations  that  relate  the  flux  densities  D,B  and  the  conductive  current  Jc  to 
the  electric  and  magnetic  fields.  We  have 

(i)  D  =eoE  +  PT2o, 

(ii)  B  =  p0H,  (2) 

,  (iii)  Jc  =  cE/fi. 


In  (2),  2fi  denotes  the  indicator  function  on  the  dielectric  medium  f 1.  Thus  Jc  =  0  in  the  ambient  or  air. 
The  total  electric  polarization  Pt  is  given  by 


Pt  =  Pi  +  P  =  cox®  +  P, 

where  Pi  is  the  instantaneous  polarization  due  to  the  interface  between  fi0  ancl  ^  and  P  is  the  material  or 
dielectric  polarization.  Hence  the  constitutive  law  (2,  i)  in  12  becomes 

D  =  eo^r-E  +  P , 


where  er  =  (1  +  %)  is  the  relative  permittivity  of  the  dielectric  medium. 

Our  main  focus  in  this  presentation  is  the  dielectric  polarization  P  which  we  assume  has  the  general 
convolution  form 


P(t,  x)  =  g  *  E(f,  x)  =  /  g(t  —  s,  x)E(s,  x)ds, 
Jo 


(3) 


where  g  is  the  general  dielectric  response  function  (DRF).  In  every  practical  example  (Debye,  Lorentz, 
etc.)  DRFs  are  parameter  dependent  as  well  as  time  (and  possibly  space)  dependent;  we  represent  this  as 
g  =  ff(2,x;  v),  where  typically  v  =  (e^ ,  es,r)  contains  parameters  such  as  the  high  frequency  limit  dielectric 
permittivity  the  static  permittivity  es,  and  relaxation  time  r.  Examples  of  often-used  DRFs  are  the 
Debye  [11,  23,  29]  in  a  material  region  fi  defined  in  the  time  domain  by 


£o(es  -  £oo )/t  e  t/r 


the  Lorentz  [11,  23,  40]  given  by 

g(t,x)  =  e0iAjl/v0e~t/2Tsin(v0t), 
and  the  Cole-Cole  [23,  27,  32,  38,  46]  defined  by 


5(f,x)  =  C  1 


Co(Cs  Cqo) 

1  + (sr)“ 


J_  fC+i°°  £0(e,-eoo) 

J Q—ioo  1  +  ( sr)a 


: Stds , 


where  C  is  the  Laplace  transform. 

These  DRFs  also  play  a  fundamental  role  in  convolution  representations  such  as  (3)  for  nonlinear  polar¬ 
ization  laws  [8,  19,  20,  24,  40]  of  Kerr  type  and  Raman  scattering  [40].  While  the  ideas  we  describe  here  on 
distributions  of  relaxation  times  and  mechanisms  can  readily  be  used  to  treat  such  nonlinear  polarization 
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laws,  we  shall  in  this  paper  concentrate  on  linear  constitutive  laws  involving  multiple  relaxation  parameters 
and  multiple  dispersive  mechanisms  in  materials  in  the  presence  of  multiple  interrogating  frequencies. 

The  macroscopic  polarization  model  (3)  can  be  derived  from  microscopic  dipole,  electron  cloud,  etc., 
formulations  by  passing  to  a  limit  over  the  molecular  population  [30].  However,  such  derivations  tacitly 
assume  that  one  has  similar  individual  (molecular,  dipole,  etc.)  parameters;  that  is,  all  dipoles,  molecules, 
“electron  clouds”,  etc.,  have  the  same  relaxation  parameters,  plasma  frequencies,  etc.  Historically,  such 
models  based  on  molecular  level  homogeneity  throughout  the  material  have  often  not  performed  well  when 
trying  to  compare  models  with  experimental  data.  Indeed,  in  1907  Von  Schweidler  [23,  50]  observed  the 
need  to  assume  multiple  relaxation  times  when  considering  experimental  data  and  in  1913  Wagner  [23,  51] 
proposed  continuous  distributions  of  relaxation  times.  This  idea  was  subsequently  visited  by  Williams  and 
Ferry  [55]  in  comparing  the  behavior  of  viscoelastic  polymers  and  dielectric  materials.  In  the  past  half 
century  intensive  experimental  efforts  [31,  32,  33,  34,  35,  38,  45]  have  been  pursued  in  describing  data  for 
complex  materials  with  distributions  of  dielectric  parameters  (especially  relaxation  times  in  multiple  Debye 
[31]  or  multiple  Lorentz  [40]  mechanisms)  in  the  frequency  domain.  A  significant  amount  of  this  work  is 
reviewed  in  the  survey  paper  by  Foster  and  Schwan  [31].  There  are  now  incontrovertible  experimentally  based 
arguments  for  distributions  of  relaxation  parameters  in  mechanisms  for  heterogeneous  materials.  Moreover, 
there  is  compelling  evidence  of  the  presence  of  multiple  mechanisms  in  complex  materials  such  as  tissue 
and  modern  polymeric  composites.  These  multiple  mechanisms  may  involve  interfacial  polarization,  dipolar 
orientation,  ionic  diffusion  (e.g.,  see  p.  40,  49,  57  of  [31])  and  may  often  require  a  selection  of  several 
types  of  distributional  representations  from  examples  such  as  the  fractional  power  laws  of  the  Cole-Cole 
[25,  26,  31,  35,  46],  the  log  normal,  the  uniform,  as  well  as  the  Debye  and  Lorentz  (although  the  fractional 
power  law  of  Cole-Cole  is  more  the  rule  rather  than  the  exception  -  p.  39,  [31]).  These  multiple  mechanisms 
are  likely  present  in  some  weighted  combination  (e.g.,  see  [36]  and  p.  369,  [40])  and  often  are  manifested  in 
a  frequency-dependent  manner.  It  is  therefore  advantageous  to  consider  interrogation  or  inverse  problems 
with  multiple  frequencies  (e.g.,  ranging  from  RF  (106)  to  GHz  (1010))  or  broadband  excitation  signals. 

To  allow  for  a  distribution  F  of  parameters  v  over  some  admissible  set  A f,  we  generalize  the  polarization 
law  (3)  to 


P(t,  x;  F)  =  h*  E(t,  x) 


o  JaT 


s,  x;  v  )E(s,  x)dF(is)ds. 


(4) 


We  expect  to  chose  F  from  (or  from  a  subspace  of)  the  space  T 


iP(A/~)  of  all  probability  measures  F  on 


a r. 


We  further  generalize  (4)  to  allow  for  dielectric  materials  with  multiple  mechanisms  or  multiple  DRFs 
(i.e.,  heterogeneous  molecular  structures)  by  considering  a  family  Q  of  possible  DRFs  and  distributions  M 
over  this  family.  This  leads  to  the  polarization  constitutive  relationship 


P(t,  x;  M,  F) 


o  Jg  JN 


—  s,x-,v)'E(s,x)dF(v)dM(g)ds  =  f  K(t  —  s,x;  M,  F)E(s,  x)ds, 

Jo 


(5) 


where  for  F  €  F  =  *P(A/")  and  M  £  M  =  K  is  defined  by 


I\{t  —  s,x;  M,  F)  =  f  f  g(t  —  s,  x;  v)dF(v)dM(g).  (6) 

Jg  J at 

When  we  use  (5)  and  (6)  in  the  Maxwell  system  (l)-(2),  we  are  led  to  a  system  of  partial  differential 
equations  where  lower  order  terms  (in  time)  depend  on  probability  measures.  These  measures  are  now  the 
“parameters”  that  characterize  the  material  dielectric  properties  which  one  must  estimate  or  identify  in 
interrogation  problems. 

With  the  recently  growing  interest  in  incorporating  uncertainty  into  models  and  systems,  the  need  to 
employ  dynamics  with  probabilistic  structures  has  received  increased  emphasis.  In  particular,  systems  with 
probability  measures  embedded  in  the  dynamics  (problems  involving  aggregate  dynamics  as  discussed  in  [8]) 
have  become  important  in  applications  in  biology  [9, 10,  8],  electromagnetics  [12]  and  hysteretic  [15,  16,  41,  42] 
and  polymeric  [17,  18,  21]  materials.  These  systems  (in  the  case  of  first  order  ordinary  differential  equations) 
have  the  form 


x(t)  =  f(t,x(t),F), 
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where  F  is  a  probability  distribution  or  measure.  In  fact  such  systems  are  not  new  and  arise  in  relaxed  or 
chattering  control  problems  [43,  44,  47,  52,  53,  54]  wherein  the  controls  are  probability  measures.  Indeed, 
such  systems  date  back  to  the  seminal  work  of  L.C.  Young  on  generalized  curves  in  the  calculus  of  variations 
[56,  57], 

In  the  next  section  we  formulate  an  inverse  problem  for  the  Maxwell  system  with  the  general  polarization 
law  (5)  and  discuss  theoretical  and  computational  aspects  of  this  problem.  We  then  report  on  our  initial 
computational  efforts  on  a  1-dimensional  version  of  the  problem  of  finding  the  underlying  polarization  in  a 
slab  of  material  using  signals  reflected  from  a  front  air/material  interface  and  a  metal  backing.  We  present 
numerical  results  for  inverse  problems  involving  Maxwell’s  equations  with  absorbing  left  boundary  condition, 
a  supraconducting  right  boundary  condition  and  a  general  polarization  term  which  includes  uncertainty  in 
the  dielectric  parameters.  We  attempt  to  determine  an  unknown  distribution  of  parameters  which  describes 
the  dielectric  properties  of  the  material.  We  explore  both  discrete  and  continuous  distributions;  for  the 
continuous  case  appropriate  parameterizations  and  discretizations  are  used. 


2  A  General  Inverse  Problem 


We  consider  the  Maxwell  system  (l)-(2)  with  polarization  P  =  P(t,  x;M,  F)  given  by  (5).  Let 

z(t  x-  M  F)  —  /^(f,  x;M,  E) 


with  (f,  x)  — >  z(t,  x;M,  F)  mapping  from  (0 ,  T)  x  SI  to  R6.  We  assume  we  are  given  data  d  =  {dj}”=1 
corresponding  to  observations  of  CJ4z(£,;,  • ;  M,  F).  Here  Ca  denotes  evaluation  of  one  or  more  components  of 
E  or  H  at  an  antenna  {x^}.  We  use  this  data  to  estimate  the  distributions  M  and  F  in  an  ordinary  least 
squares  (OLS)  formulation,  seeking  to  minimize 


J(M,  F)  =  Y,  I Ca^U,  ■ ;  M,  F)  -  d8|2  (7) 

2  =  1 


over  (. M,F )  G  M.  x  T  =  ty(G)  x  fP(AT)-  We  thus  seek  to  find  (. M,F )  such  that 

(M,  F)  =  arg  min {J(M,  F)  :  M  e  M,  F  €  T}. 


We  note  that  while  for  simplicity  we  use  an  OLS  formulation  here,  most  of  the  results  discussed  below 
could  readily  be  developed  in  the  context  of  other  standard  estimation  formulations  such  as  maximum 
likelihood  estimators  (MLE),  weighted  least  squares  (WLS),  or  generalized  least  squares  (GLS). 

For  theoretical  and  computational  purposes,  one  needs  a  topology  on  A4  and  T  and  for  this  we  choose  the 
Prohorov  metric  p*  of  weak*  convergence  in  A4  and  T  when  they  are  considered  as  subsets  of  the  topological 
duals  C%{Q)  and  C*B (J\f)  of  the  spaces  Cb(G)  and  C'b(AT)  of  bounded  continuous  functions  on  Q  and  A/", 
respectively  [5,  8].  That  is,  Fk  — >  F  in  the  p*  metric  if  and  only  if 


(j>{v)dFk{v) 


<j>(v)dF(v) 


for  all  (f>  €  C'b (AT) ,  i.e.,  all  bounded  continuous  <f>  on  A/";  similarly  for  Mk  — >  M.  It  is  known  [5,  8]  that  if 
G  and  A f  are  complete  metric  spaces,  then  A4  =  ty(G)  and  T  =  ^(A/)  taken  with  the  Prohorov  metric  p* 
are  complete  metric  spaces.  Moreover,  if  Q  and/or  Af  are  compact,  then  so  are  A4  and/or  T .  Using  these 
properties  and  arguments  similar  to  those  in  [5,  14],  the  following  problem  stability  results  can  be  proven. 


Theorem  1:  Suppose  ( M,F )  — *  Ca z(t,  •  ;M,F)  is  continuous  on  M.  x  T  and  suppose  that  Q  and  M 
are  compact.  Then  solutions  ( M,F )  of  minimizing  (7)  exist  (generally,  non-uniquely)  and  are  continuous  in 
the  data  d  in  the  following  sense.  Suppose  (M*(d), F*(d))  and  (M*(dfc), F*(dk))  are  the  solution  sets  of 
minimizing  (7)  for  data  d  and  dfc,  respectively,  where  dfc  — >  d  as  k  — >  oo.  Then 

dist  [(M*(dk),F*(dk))  ,  (M*(d),  F*(d))]  -►  0 
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as  k  — >  oo,  where  dist[  ,  •  ]  is  the  Hausdorf  distance  in  the  metric  space  M  x  T. 

The  problem  of  minimizing  (7)  over  M  x  T  is  in  general  an  infinite  dimensional  optimization  problem 
that  poses  formidable  computational  challenges.  But  we  note  that  set  Q  of  possible  DRFs  can  in  almost  all 
cases  be  taken  as  a  fixed  finite  set,  i.e.,  Q  =  {<71,32,  •••,  <?A'}  which  means  that  M.  would  be  given  by 

M  =  \  M(g)  G  <£(£)  |  dM{g)  =  ^  cy  dAgj  ( g),aj  >  0,  ^  aj  =  1 
[  3= 1  3 

where  Ag.  is  the  Dirac  measure  with  single  atom  at  gj,  i.e.,  dAg  (g)  =  5gj(g)dg.  Thus,  this  part  of  the 
optimization  reduces  to  one  over  a  closed  bounded  convex  set  in  Euclidean  space.  The  minimization  over 
T  =  ip [M )  is  more  interesting  since  in  general  one  expects  A f  to  involve  a  continuum  of  (vector)  parameters. 
We  illustrate  the  possibilities  with  Af  =  T  =  {  r  |r  G  [ra,T(,]}  for  r  the  relaxation  parameter  in,  for  example, 
a  Debye  or  Lorentz  mechanism.  There  are  a  number  of  ways  this  problem  could  be  approached: 

1.  Assume  that  F  is  discrete,  having  the  fixed  form  dF{r)  =  J2ajST At) dr  and  seek  to  find  ( aj,Tj ) 
minimizing  (7)  where  ay  >  0 =  1  ,  D  ^  [ TaiTb\\ 

2.  Assume  that  F  is  (absolutely)  continuous  and  given  in  parametric  form  dF{r)  =  f{T,p,a)dr 
where  f  is  known  (e.  g.,  normal,  log-normal,  uniform,  etc.)  and  seek  to  find  (/ u,a ); 

3.  Assume  F  does  not  have  a  specific  parametric  form  and  seek  to  find  the  general  form  for  F  through 
the  optimization  of  (7). 

In  the  first  two  cases  above,  one  effectively  reduces  the  inverse  problem  to  a  computationally  tractable 
(one  hopes  !-we  explore  these  ideas  computationally  in  subsequent  discussions  below)  optimization  problem 
that  is  finite  dimensional.  The  third  case  remains  infinite  dimensional  in  nature  and  one  must  develop 
approximation  ideas  that  lead  to  implementable  computational  algorithms.  In  [5],  the  authors  developed 
approximation  ideas  based  on  density  results  for  measures  arising  in  probability  theory.  We  only  outline 
those  here,  referring  the  reader  to  [5]  for  more  details  and  proofs. 

To  develop  approximation  ideas  for  the  nonparametric  case,  we  first  consider  a  family  of  partition  points 
TN  =  { qj N  =  1,2, ...,  such  that  U^’=1TAr  is  dense  in  T.  Then  define 

N 

jrN  =  spWj'T')  =  {pN  e  f$(T)\dFN (t)  =  ^(r)dr,  qf  G  TN,pf  rational,  =  1}. 

3= 1 

It  can  be  argued  that  U^,_1lFAr  is  dense  in  T  in  the  Prohorov  metric  p*.  Moreover,  if  T  is  compact, 
one  can  prove  a  method  stability  theorem  (see  [5,  8])  similar  to  Theorem  1  above.  Specifically,  let  F^( d)  be 
the  set  of  solutions  obtained  in  minimizing  (7)  with  T  replaced  by  TN .  Then  the  method  stability  theorem 
guarantees 

dist[(M*(dk),F^(dk)),(M*(d),F*{d))}  -  0 

as  N,  k  — >  00. 

More  generally,  one  may  wish  to  consider  only  classes  of  distributions  T  that  arise  from  densities  functions, 
i.e.,  absolutely  continuous  distributions  in 

Tag  =  {  F  G  T  \  F'  =  /,  /  G  Fweak} 

where  FU!eoj-  is  a  given  weakly  compact  subset  of  L2(T).  For  example,  Fweak  could  be  any  given  closed, 
convex,  bounded  subset  of  L2(T).  It  is  proven  in  [21]  that  sets  T ac  defined  in  this  way  are  compact  (in 
the  p*  metric)  subsets  of  T  and  the  corresponding  existence  and  stability  results  of  Theorem  1  hold  for 
problems  using  Tag  in  place  of  T  in  the  optimization  for  (7).  But  once  again  these  are  infinite  dimensional 
in  nature  and  approximations  are  needed.  We  note,  that  although  the  computational  framework  described 
above  utilizing  Dirac  measures  is  also  valid  here,  it  is  often  desirable  to  develop  “smoother”  approximations 
to  elements  of  T ac-  In  particular,  suppose  that  /  G  Tweak  and  F  G  T  =  %TT)  with  F  =  f  f.  Since  Tweak 
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is  a  subset  of  L2(T),  we  can  use  any  number  of  types  of  splines  to  formulate  an  approximation  to  /,  i.e.,  let 
{Sj'jjLi,  N  =  1,2,...  be  a  standard  spline  family  [14,  48,  49]  (e.g.,  piecewise-linear,  cubic,  Hermite  cubic, 
etc.)  and 

N 

/"M  =  £<sf(  r), 

3= 1 

where  the  rational  numbers  uA  are  chosen  so  that  fN  — >  /  in  L2(T).  This  implies  that 


(t>fdr 


for  all  <j)  £  L2(T ),  and  hence  for  all  4>  £  C(T ),  which  yields 


p(FN ,  F)  — >  0, 


where  FN  =  f  fN .  Thus  defining 


N 


?SPL  ={h£  L2{T)\h{r)  =  Y,wiSi(r)  )  ’ 

f=i 


we  can  conclude  that  the  set 


Tag  =  {F£F\F  =  J  fj  £  U%=1FgPL} 

is  dense  in  Fa c  in  the  p* -metric.  Hence  spline  families  provide  an  alternative  way  to  approximate  elements 
of  F ac  in  our  computational  work.  We  note  that  in  this  case  the  approximating  elements  are  not  probability 
distributions  themselves.  Once  again  one  can  state  and  prove  a  method  stability  theorem  [14]  using  these 
spline  approximations.  That  is: 

Theorem  2:  Suppose  T  is  compact  and  the  solutions  C^z(i,  •  ;M,F )  are  continuous.  Let  Fweak,  F ac , 
and  FgPL  be  as  defined  above  with  (M*(d),  F^(d))  the  solution  sets  for  minimizing  (7)  with  F  replaced  by 
FgpL •  Then  for  dfc  — >  d  we  have 

dist[(M*(dk),F^(dk)),(M*(d),F*(d))]  ->  0 

as  N,  k  — »  oo.  Thus  solutions  depend  continuously  on  data  and  the  approximate  problems  are  method  stable. 


To  illustrate  computational  aspects  of  the  ideas  presented  in  this  section,  we  turn  next  to  a  1-D  version 
of  the  inverse  problem  and  present  results  in  the  next  several  sections  on  use  of  discrete  and  continuous 
distributions  in  the  polarization  laws. 


3  The  1-D  Problem  Formulation 

For  our  initial  numerical  efforts,  we  turned  to  the  1-D  example  as  explained  in  detail  in  [11].  Under  the 
assumptions  detailed  there,  one  obtains  a  domain  as  depicted  in  Figure  1. 

Restricting  to  one  dimension,  and  using  D  =  eE  +  P,  we  can  write  Maxwell’s  equations  in  second  order 
form  as 

p,0eE  +  poInP  +  P-o crdnE  -  E"  =  —p,0  js  in  fi  U  fi0,  (8) 

where  E  is  the  k  or  2  component  of  the  electric  field,  P  is  the  media’s  macroscopic  electric  polarization,  Js 
is  the  interrogating  signal,  f 1  is  the  domain  of  the  material  under  investigation,  fio  is  the  ambient  domain 
(considered  a  vacuum),  and  a  =  cr(z)  is  the  conductivity  of  the  material.  Note  we  are  considering  a  non¬ 
magnetic  material  containing  no  charge  distribution  (p  =  0).  Also,  let  e  =  eo(l  +  Ln(er  —  1))  where  er  =  er(z) 
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z  (meters) 


Figure  1:  The  1-D  domain. 


is  the  dielectric  permittivity.  The  values  /ro  and  eo  are  the  magnetic  permeability  and  the  electric  permittivity 
of  free  space,  respectively.  See  [11]  for  more  details.  The  boundary  conditions  that  we  are  assuming  are 
absorbing  at  z  =  0  (to  provide  a  finite  computational  window)  and  supraconducting  at  z  =  1  (representing 
a  metal  backing) : 


E-cE' 


=  0 

J  z= 0 

E(t,l)  =  0. 


We  assume  homogeneous  initial  conditions 


E{0,z)  =  0 
E(0,z)  =  0. 


(9) 


(10) 

(11) 


For  implementation,  we  scale  time  by  t  =  ct  and  the  polarization  by  P  =  P/e o  for  convenience.  Also 
note  that  we  have  employed  the  “method  of  mappings”  to  obtain  a  computational  domain  oi  z  =  [0,1].  The 
actual  dimensions  of  the  domains  considered  in  this  report  depend  on  the  interrogating  frequency  u>.  In 
particular,  for  u>  =  2n  x  1011  we  consider  a  material  slab  of  thickness  .004?n  preceded  by  a  vacuum  of  depth 
.002m.  We  scale  these  dimensions  as  the  wavelength  scales  when  we  change  interrogating  frequencies. 
Converting  the  system  (8)  to  weak  form  (and  dropping  the '  notation)  we  obtain 

(erE(t,  -),4>)  +  (InP(t,  -),4>)  +  {voaInE(t,  •),</>)  -  (E"(t,  •),</>)  =  -(rj0js(t,  •),</>), 

where  rj o  =  \J /Jq / £o  and  </>  G  V  =  ^^(0, 1)  =  {4>  £  H1( 0, 1)  :  ^>(1)  =  0}.  Finally,  we  integrate  by  parts,  and 
apply  the  boundary  conditions  (9)  to  get 

(erE(t,  ■),</))  +  (InP(t,  •),</>)  +  (r)o(rInE(t,  ■),</))  +  ( E'(t ,  •),</>')  +  E(t,  0)^(0)  =  -( rj0js(t ,  •),</>)•  (12) 


To  describe  the  behavior  of  the  electric  polarization  P,  we  may  employ  a  general  polarization  kernel,  or 
dielectric  response  function,  g  as  follows: 


P(t,z)=  /  g(t  —  s,  z\t)E(s,  z)ds 

J  o 


(13) 


7 


where,  for  instance  using  a  Debye  polarization  model, 


9(t\  t) 


eo(es  -  eoo )/t  e  t/r . 


However,  this  presupposes  that  the  material  may  be  sufficiently  defined  by  a  single  relaxation  parameter 
r,  which  is  generally  not  the  case.  In  order  to  account  for  uncertainty  in  the  polarization  mechanisms,  we 
allow  for  a  distribution  of  relaxation  parameters.  Thus,  we  define  our  polarization  model  in  terms  of  a 
distribution  dependent  dielectric  response  function  h 


P(t,z)=  (  h(t  —  s,  z\  F)E(s,  z)ds, 

J  o 


(14) 


where  h  is  determined  by  a  family  of  polarization  laws  each  described  by  a  different  parameter  r,  and 
therefore  is  given  by 

h(t,z;F)  =  J  g(t,  z;  r)dF(r), 

where  T  =  [ti,  t2]  C  (0,  oo).  In  particular,  if  the  distribution  F  were  discrete,  consisting  of  a  single  relaxation 
parameter,  then  we  would  again  have  (13). 

The  macroscopic  electric  polarization  becomes 


P{t,z) 


—  s,  z\  r)dF{r) 


E(s,  z)ds , 


or,  interchanging  integrals,  we  have 


p(tiz)  =  J  'P(t,z-,T)dF(r), 


where 


V(t,z\r)=  [  g(t  —  5,  z\  t)E(s,  z)ds 
Jo 


is  the  polarization  due  to  the  relaxation  parameter  r.  Thus  assuming  we  have  a  computational  method  to 
compute  (13),  i.e.,  V,  we  use  this  as  a  basis  for  approximating  P,  either  directly  in  the  discrete  case,  or  using 
a  quadrature  rule  in  the  continuous  case. 

The  theory  developed  in  Section  2  for  general  inverse  problems  can  be  directly  applied  to  this  1-D  problem 
once  one  has  established  continuity  of  the  solution  E  with  respect  to  the  measures  F  in  the  space  T  =  *P(T) 
taken  with  the  Prohorov  metric.  But  these  desired  continuity  results  are  given  in  [12] . 


3.1  Discrete  Distribution 

Consider  the  discrete  distribution  given  by 


dF(r) 


E 


d(n) 

£+1 


dr , 


where  S  represents  the  Dirac  distribution,  and  Tt  =  ra  +  ir h  with  =  Tb  j Ta  .  Then  we  have 


e 

P(t,z)  =  ’Y^aiV(t,z\Ti), 

i=0 

where  a*  =  C  j .  Note  that  for  this  example,  the  nodes  (t;)  are  linearly  spaced  and  the  weights  or  masses 
(ai)  are  uniform.  More  generally,  one  can  treat  a  discrete  distribution  without  linearly  spaced  nodes  and/or 
uniform  masses. 


3.2  Uniform  Distribution 

The  simplest  continuous  distribution  is  a  uniform  distribution.  Consider 

dF{r )  =  - dr, 

T~b  T~a 

for  ra  <  t  <  Tb ,  and  zero  otherwise.  Since  the  distribution  function  is  constant,  we  may  pull  it  outside  the 
integral  giving 

P(t,z)= — — —  [  V(t,z;r)dT.  (15) 

Z}f  Ta  J  Ta 

While  this  reduces  down  to  essentially  integrating  r  out  of  the  polarization  convolution  term,  for  most  polar¬ 
ization  models  this  does  not  have  an  analytical  solution.  Therefore  we  must  resort  to  numerical  quadrature. 
Recall  the  Composite  Simpson’s  rule  approximation  to 

f [«,&](£)  :=  [  L(x)dx 


IS 


2=0 


where  l  is  even,  Xi  =  a  +  ih ,  h  =  — ,  and  the  weights  are  given  by 


(  |  if  *  =  0  or  *  = 


c,  =  <  |  else  if  i  odd 


|  else  if  i  even 


Thus,  our  Composite  Simpson’s  approximation  to  (15)  can  be  written 


1  £ 

P(t,  z)  =  - V  4V(t,  z ;  Ti)Th 

Tb  ~  T“  U 

l 

~  OLjPjt,  z\  Tj) 


2=0 


with 


since  the  numerator  in  t>,  =  Tb^Ta  cancels  the  constant  in  front  of  the  integral.  Note  here  that  while  our 
nodes  are  still  linearly  spaced,  now  our  weights  are  non-uniform.  We  could  have  used  a  uniform  discrete 
distribution  to  approximate  the  continuous  one,  which  would  have  resulted  in  uniform  weights,  but  the 
corresponding  quadrature  rule  for  g(t,  z:  r)  would  have  only  been  0(jh)-  Composite  Simpson’s  rule  provides 
better  accuracy  for  the  same  number  of  function  evaluations. 


3.3  Log-Normal 

A  more  realistic  model  for  the  distribution  of  Debye  relaxation  times  is  given  by  the  log-normal  distribution 
(see  [23]).  Therefore  we  consider  the  probability  distribution  defined  by 


dF(r) 


1  1 


ln  10  T 


exp 


(log  t  /r) 2  \ 
2  cr2  J 


dr, 
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which  means  that  log  r  is  normally  distributed  with  mean  fj,  and  variance  a2.  The  In  10  term  appears  because 
we  are  using  base  10  logarithms.  Thus 


P(t,z)=  /  V(t,  z\  t)cIF(t) 


t)oLF(t), 


(16) 


for  r a  sufficiently  close  to  zero  and  tj,  sufficiently  large.  Since  logr  is  normally  distributed,  we  choose  our 
nodes  based  on  a  uniform  discretization  of  logr,  denoted  {£i}i=o  from  logr0  :=  /j  —  6cr  to  log rt,  :=  /j  +  6tr, 
where  ra  and  r i,  are  thusly  defined  by  their  logarithms.  In  this  way  we  reduce  our  support  to  a  compact  set 
which  should  be  large  enough  to  effectively  approximate  our  integral.  In  order  to  use  Composite  Simpson’s 
rule  we  prefer  to  have  a  uniform  discretization  of  nodes.  Therefore  we  first  change  variables  so  that  the 
integral  is  in  terms  of  £  =  log  r.  First,  let 


/(o  =  ^W  (logT^ 


7T<7^ 


2a2 


and  note  that 


Then  we  have  that 


d(  =  d 


1  -dr. 


P(t,  z) 


In  r 

In  10  )  In  10  t 


P(t,z ;  r)dF(r)  = 


'  fi—6a 


P(t,z ;  10  «)/m 


where  implicitly  dF  =  f{r)dr  and  dF  =  f(£)d£.  Applying  Composite  Simpson’s  rule  gives  finally 


i 

P(t,  z)  «  ^2  Z\  Ti) 

i= 0 


where  we  define  t,  =  10^  and 


ai  =  4f(Q£h  = 


cf  12cr  1 


(■  V2 7 


:  exp  -- 


(logn  - 


2a2 


For  this  case  our  nodes  are  non-linearly  spaced  in  t  and  our  weights  are  non-uniform. 

The  method  used  in  the  bi-gaussian  case  is  derived  in  a  similar  fashion.  Consider  that  the  polarization  is 
driven  by  two  distinct  mechanisms,  one  with  a  dielectric  distribution  determined  by  mean  fj>i  and  standard 
deviation  <j\,  and  the  other  determined  by  mean  /i2  and  standard  deviation  cf2-  Then  we  define  F\  (r;  /xi,cri) 
and  F2(t;  fi 2,(72)  as  log  normal  distributions  as  above,  and  let  our  macroscopic  electric  polarization  be  a 
function  of  some  combination  of  these  distributions  (e.g.,  determined  by  the  relative  volume  percentage 
of  each  of  two  substances  in  a  material).  Thus  if  dF(r)  =  /3idFi(r)  +  (hdF-iiT)  then  we  again  have  the 
representation  (16).  For  the  discretization,  we  prefer  to  apply  Composite  Simpson’s  to  each  distribution 
separately  and  then  combine  at  the  end.  In  other  words 

i  i 

P(t,  z)  «  /3 a}V(t ,  r/)  +  /32  ^  ^  Ti )> 

i— 0  i—O 

where  {a}}  and  {r? }  are  determined  by  ft\  and  <71,  and  {a2}  and  {t2}  are  determined  by  p2  and  c72,  as 
explained  above. 


4  Inverse  Problem  Formulation 

Our  goal  is  to  estimate  the  probability  distribution  function  (PDF)  F  £  'P(T)  of  relaxation  parameters 
in  a  given  model  of  the  polarization,  where  V(T)  is  the  set  of  all  PDFs  on  the  admissible  region  T  = 
[t-! , 7”2 ]  C  (0,oo),  by  using  reflections  of  electromagnetic  interrogating  signals  off  a  metallic  backing  of  a 
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dielectric  material.  To  this  end  we  attempt  to  minimize  the  difference  between  simulations  and  observations 
of  time-domain  data.  In  our  formulation,  the  observations,  Ej,  are  of  the  electric  field  E  at  discrete  times 
tj  taken  at  z  =  0.  Each  simulation  is  a  solution  of  Maxwell’s  equations  given  in  (12)  using  the  polarization 
model  (14)  with  candidate  values  for  the  distribution  of  relaxation  parameters.  We  propose  a  non-standard 
least-squares  measurement  (see  [13]  for  a  detailed  justification)  for  the  objective  function,  which  is  given  by 

W  =  E|  (17) 

3 

where  E( •,  •;  F)  is  the  solution  to  (12)  and  (14)  corresponding  to  the  distribution  F.  Thus  the  inverse  problem 
is  to  solve 

min  J(F). 

FGV(T) 

5  Inverse  Problem  Results  Using  a  Discrete  Distribution 

Consider  a  discrete  distribution  with  two  atoms  at  T\  and  T2-  Essentially,  we  are  decomposing  the  electric 
polarization  into  two  components,  each  dependent  on  distinct  relaxation  times  as  follows: 

P  =  a{P{t,  z;n)  +  a2V(t,z-,T2),  (18) 

where  each  satisfies  a  Debye  polarization  equation  with  parameter  iy.  For  now  we  assume  the 

proportions  aq  and  =  1  —  cq  are  known.  Thus  we  are  attempting  to  solve  the  following  least  squares 
optimization  problem: 

,mil\5Zl  IUU0;('ri,T2))|  -  \Ej\  ,  (19) 

(r i  ,t2  )  I 

3 

where  Ej  is  synthetic  data  generated  using  the  true  solution  (r* ,  )  in  our  simulator  (with  a  highly  refined 

mesh),  and  E(tj,  0;  (71,72))  depends  on  each  r,;  through  its  dependence  on  P,  see  for  example  (12).  Figures 
2  and  3  depict  an  example  of  the  objective  function  and  the  log  of  the  objective  function  respectively,  both 
plotted  versus  the  logs  of  n  and  T2  (using  a  frequency  of  1011  Hz,  oq  =  a  2  =  -5,  rj*  =  ICC7,5  and  r|  =  10  7‘8) . 

5.1  Analysis  of  Objective  Function 

We  can  see  clearly  from  the  log  surface  plot  in  Figure  3  that  there  exists  a  relation  for  the  relaxation  times 
for  which  the  corresponding  simulations  best  match  the  data.  We  will  refer  to  this  relation  as  the  “curve  of 
best  fit”. 

Note  that  the  appearance  of  many  local  minima  is  due  to  the  steep  decent  near  the  “curve  of  best  fit” 
since  the  lattice  points  of  the  mesh  used  do  not  always  lie  near  the  curve.  If  we  trace  along  this  curve,  as 
displayed  in  Figure  4,  we  see  that  there  are  actually  two  global  minima,  the  exact  solution  of  log(ri)  =  —7.5 
and  log(r2)  =  —7.8,  and  since  the  proportion  used  in  this  case  was  oq  =  0:2  =  -5,  we  also  have  the  symmetric 
solution  where  Ti  and  72  are  swapped.  If  we  superimpose  the  data  from  the  contour  plot  and  the  data  from 
the  surface  plot  on  a  lattice,  as  demonstrated  in  Figure  5,  we  can  more  easily  see  the  structure  of  the  “curve 
of  best  fit” ,  and  both  global  minima. 

Unfortunately,  the  scale  of  these  plots  shows  that  the  difference  between  the  objective  function  at  the 
minimizers  and  any  other  point  on  this  curve  in  our  parameter  space  is  less  than  4  x  10-10.  Therefore  the 
exact  minimizing  parameters  are  not  likely  to  be  identifiable  in  a  practical,  experimental  setting. 

The  equation  for  this  “best  fit  curve”  can  be  derived  by  combining  equations  (12)  and  (14).  Note  that 

P{t,z)  =  f  Q(t  —  s,  z)E(s,z)ds  +  Q(0,z)E(t,z)  +  Q(0,z)E(t,z),  (20) 

J  0 

where 

G(t,z;F)  =  J  g(t,  z\t)cLF{t). 

So  that  now  substituting  (20)  into  (12)  we  obtain 
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f=1e11,a  =.5 


Figure  2:  The  objective  function  for  the  relaxation  time  inverse  problem  versus  the  log  of  r\  and  the  log  of 
r2  using  a  frequency  of  10 nHz. 

f=1e11,a1=.5 


Figure  3:  The  log  of  the  objective  function  for  the  relaxation  time  inverse  problem  versus  the  log  of  T\  and 
the  log  of  r2  using  a  frequency  of  lQxlHz.  The  solid  line  above  the  surface  represents  the  curve  of  constant 
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f=1  el  1,lt=. 200556 


Figure  4:  The  objective  function  for  the  relaxation  time  inverse  problem  plotted  along  the  curve  of  constant 
A  using  a  frequency  of  IQ11  Hz. 


Suface  plot  superimposed  with  contour  plot 


Figure  5:  The  surface  plot  of  the  objective  function  for  the  relaxation  time  inverse  problem  superimposed 
with  the  contour  plot  along  the  curve  of  constant  A,  using  a  frequency  of  10nHz.  The  fact  that  there  are 
only  two  minima  is  clearly  more  visible  in  this  plot. 
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{erE(t,  ■),(/))  +  {In  [vo<r  +  G{ 0,  z)\  E(t,  •),</>)  +  (InG( 0,  -)E(t,  •),<£) 

+  (  [  InG(t  -  s,  -)E(s,  ■ )ds ,  <t>)  +  {E'(t,  •),</>')  +  £(i,  0)0(0)  =  -(r)0js(t,  •),  0)-  (21) 

Jo 

For  large  frequencies,  the  E  term  in  (21)  dominates  over  all  other  terms  which  depend  on  Q .  We  can  see 
this  by  considering  the  frequency  domain  where  the  time  derivative  causes  an  increase  by  a  factor  of  u>  in 
the  norm.  Thus  the  equation  for  the  “curve  of  best  fit”  is  simply  that  of  constant  1/(0,  z).  For  the  discrete 
Debye  example  we  have 


G(0,z)  =  J  ff(0,2;r)dF(r)  =  aie0(es  -  £oo)/n  +  a2e0(es  -  eoo)/T2, 
therefore  the  equation  is 


Of  1  Q?2  OL 1  OL  2 


Tl  T2 


'1  r2 

This  is  precisely  the  curve  that  is  plotted  above  the  surface  in  Figure  3. 

For  scaling  purposes,  and  to  be  consistent  with  the  presentation  in  [11],  we  define 


A  := 


1 


^o(^s  ^oo) 


S(0,z), 


(22) 


which  in  this  example  means 

^  _  op  «2 
cri  cr2 

and  thus  we  may  say  that  the  “curve  of  best  fit”  is  the  line  of  constant  A. 

The  frequency  dependence  of  the  E  term  in  (21)  suggests  that  for  smaller  frequencies  it  may  not  be  the 
dominant  contributor,  and  therefore,  there  may  be  a  fundamentally  different  structure  to  the  surface  plot. 
This  is  in  fact  what  we  observe  in  our  simulations.  Figures  6  through  9  display  the  surface  plots  and  the  log 
surface  plots  for  frequencies  of  10 9 Hz  and  10 6 Hz.  Note  that  in  the  latter  case  the  concavity  of  the  “curve 
of  best  fit”  has  swapped  orientation!  This  demonstrates  that  the  surface  plots  are  very  much  dependent  on 
ui  even  though  the  relaxation  mechanisms  are  the  same. 

Through  our  numerical  calculations  we  have  determined  that  for  the  case  using  a  frequency  of  106  Hz  the 
“curve  of  best  fit”  is  actually  that  of  constant  t  :=  oq t\  +  a2r2  which  is  what  one  might  expect  as  this  is  the 
weighted  average  of  the  relaxation  times.  For  the  example  given  here,  f  ss  2.37000  x  10-8.  In  our  simulations, 
angular  frequencies  less  than  one  divided  by  this  number  were  characterized  by  a  constant  r  while  larger 
frequencies  resulted  in  dominance  by  the  A  term  given  above.  The  fact  that  the  regime  characterized  by 
lot  <  1  is  fundamentally  different  in  many  respects  from  that  of  the  lot  >  1  regime  is  well  documented  (see, 
for  example,  [23]).  Still,  the  behavior  of  the  objective  function  along  its  corresponding  “curve  of  best  fit”  is 
similar  for  each  frequency,  despite  the  curves  themselves  being  fundamentally  different,  as  demonstrated  by 
comparing  Figure  4  to  Figure  10  which  uses  a  frequency  of  10 6 Hz.  Note,  however,  that  the  scale  of  Figure 
10  is  several  orders  of  magnitudes  larger,  suggesting  that  the  global  minimizers  may  be  easier  to  find  for 
smaller  frequencies. 

Remark  1  We  should  mention  here  that  in  order  to  analyze  only  the  relationship  between  the  interrogating 
frequency  and  the  relaxation  time,  we  scaled  the  dimensions  of  the  interrogated  object  according  to  the  change 
in  scale  of  frequency.  For  example,  when  the  frequency  is  divided  by  100  to  go  from  1011  to  10  ,  we  accordingly 
multiply  the  slab  thickness  by  100  to  get  Am  in  order  that  the  ratio  of  the  dimensions  to  the  wavelength  of 
the  signal  remains  the  same. 

5.2  Optimization  Procedure  and  Results 

We  attempt  to  apply  a  two  parameter  Levenberg-Marquardt  optimization  routine  to  the  modified  least 
squares  error  between  the  given  data  and  our  simulations,  as  defined  in  (17),  to  try  to  identify  the  two  distinct 
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f=1e9,a1=.5 


Figure  6:  The  objective  function  for  the  relaxation  time  inverse  problem  versus  the  log  of  t\  and  the  log  of 
r2  using  a  frequency  of  10 9 Hz. 


f=1e9,a1=.5 


Figure  7:  The  log  of  the  objective  function  for  the  relaxation  time  inverse  problem  versus  the  log  of  T\  and 
the  log  of  r2  using  a  frequency  of  10 9 Hz.  The  solid  line  above  the  surface  represents  the  curve  of  constant 


A. 
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f=1e6,a1=.5 


Figure  8:  The  objective  function  for  the  relaxation  time  inverse  problem  versus  the  log  of  T\  and  the  log  of 
r2  using  a  frequency  of  10 6 Hz. 

f=1e6,a  =.5 


Figure  9:  The  log  of  the  objective  function  for  the  relaxation  time  inverse  problem  versus  the  log  of  T\  and 
the  log  of  72  using  a  frequency  of  106 Hz.  The  solid  line  above  the  surface  represents  the  curve  of  constant  t . 
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f=1e6,tt=2.37E-8 


Figure  10:  The  objective  function  for  the  relaxation  time  inverse  problem  plotted  along  the  curve  of  constant 
t  using  a  frequency  of  10 6 Hz. 


relaxation  times  that  generated  the  data.  We  are  assuming  that  the  corresponding  volume  proportions  of 
the  two  materials  (oq  and  a2  :=  1  —  aq)  are  known.  We  consider  three  different  scenarios  with  respect  to  the 
volume  proportions:  aq  £  { .  1 ,  .5,  .9} .  We  also  perform  our  inverse  problem  using  the  frequencies  1011Hz, 
109Hz  and  10 GHz.  Lastly  we  test  the  optimization  procedure  with  three  various  initial  conditions  given  in 
Table  1.  The  actual  values  are  n  =  icr7-50031  «  3.16  x  10"8  and  r2  =  icr7-80134  w  1.58  x  10"8.  Note  that 
the  first  set  of  initial  conditions  is  the  farthest  from  the  exact  solution,  while  the  third  is  the  closest. 


Table  1:  Three  sets  of  initial  conditions  for  the  relaxation  time  inverse  problem  representing  = 

(Ct^ /C)  for  C  £  {5,2,1.25}  respectively  (case  0  represents  exact  solution),  also  given  are  the  log10  of 
each  relaxation  time,  as  well  as  the  %  relative  error  from  the  exact  value. 


case 

Tl 

t2 

log(Ti) 

log(r2) 

%  Tl 

%  t2 

0 

3.1600e-8 

1.5800e-8 

-7.50031 

-7.80134 

0 

0 

1 

1.5800e-7 

3.1600e-9 

-6.80134 

-8.50031 

400 

80 

2 

6.3200e-8 

7.9000e-9 

-7.19928 

-8.10237 

100 

50 

3 

3.9500e-8 

1.2640e-8 

-7.40340 

-7.89825 

25 

20 

The  results  of  the  optimization  are  given  in  Tables  2  through  7.  Most  of  the  cases  appear  not  to  have 
converged.  Only  a  few  of  the  cases  corresponding  to  the  closest  initial  conditions  converged  close  to  the 
original  values  of  the  relaxation  times.  But  if  we  recall  the  shape  of  the  objective  function  that  we  are  trying 
to  minimize  (for  example  see  Figures  2  and  3)  then  it  is  understandable  how  a  gradient  based  method  may 
converge  directly  to  the  “curve  of  best  fit”.  However,  once  it  reaches  this  curve,  the  optimization  routine 
may  jump  back  and  forth  across  the  “trench”  and  may  not  be  able  to  traverse  the  curve  to  find  the  true 
global  minimum.  This  is  in  fact  what  is  occurring. 
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Table  2:  Resulting  values  of  T\  from  the  Levenberg-Marquardt  routine  using  a  frequency  of  10nf/z  (recall 
the  exact  solution  r*  =3.1600e-8). 


case 

.1 

ai 

.5 

.9 

i 

1.57218e-07 

1.7488e-07 

2.26538e-07 

2 

6.34973e-08 

9.73445e-08 

4.90169e-08 

3 

3.98062e-08 

4.13321e-08 

3.43569e-08 

Table  3:  Resulting  values  of  T\  from  the  Levenberg-Marquardt  routine  using  a  frequency  of  109ii2  (recall 
the  exact  solution  =3.1600e-8). 


case 

.1 

a  i 

.5 

.9 

i 

1.60128e-07 

1.60397e-07 

7.0395e-08 

2 

6.37815e-08 

6.68561e-08 

3.99819e-08 

3 

3.97416e-08 

4.03239e-08 

3.59503e-08 

Table  4:  Resulting  values  of  iq  from  the  Levenberg-Marquardt  routine  using  a  frequency  of  10 eHz  (recall 
the  exact  solution  t{  =3.1600e-8). 


case 

.1 

.5 

.9 

i 

3.79957e-08 

3.23393e-08 

3.16236e-08 

2 

3.17753e-08 

3.32218e-08 

3.21036e-08 

3 

3.17897e-08 

3.19001e-08 

3.16068e-08 

Table  5:  Resulting  values  of  T2  from  the  Levenberg-Marquardt  routine  using  a  frequency  of  lQllHz  (recall 
the  exact  solution  r|  =1.5800e-8). 


case 

.1 

Ot,\ 

.5 

.9 

i 

1.51271e-08 

1.1208e-08 

3.2429e-09 

2 

1.53697e-08 

1.18119e-08 

6.10283e-09 

3 

1.56197e-08 

1.41322e-08 

1.1607e-08 

Table  6:  Resulting  values  of  T2  from  the  Levenberg-Marquardt  routine  using  a  frequency  of  10 9Hz  (recall 
the  exact  solution  =1.5800e-8). 


case 

.1 

0^1 

.5 

.9 

i 

1.51255e-08 

1.12739e-08 

4.54532e-09 

2 

1.53692e-08 

1.25029e-08 

8.12599e-09 

3 

1.56222e-08 

1.4257e-08 

1.02271e-08 

Table  7:  Resulting  values  of  T2  from  the  Levenberg-Marquardt  routine  using  a  frequency  of  10 6Hz  (recall 
the  exact  solution  r|  =1.5800e-8). 


case 

.1 

.5 

.9 

i 

1.51064e-08 

1.50634e-08 

1.63437e-08 

2 

1.57826e-08 

1.41845e-08 

1.12902e-08 

3 

1.57792e-08 

1.55032e-08 

1.57796e-08 
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Table  8  displays  the  values  of  A  corresponding  to  each  set  of  initial  conditions,  where  case  0  denotes  the 
exact  solution.  Tables  9  and  10  display  the  resulting  values  of  A  from  the  Levenberg-Marquardt  routine 
using  frequencies  lO^Hz  and  10 9Hz  respectively.  Since  smaller  frequencies  are  expected  to  converge  to  the 
curve  of  constant  f,  we  give  the  initial  values  of  f  as  well,  in  Table  11,  and  Table  12  displays  the  resulting 
values  of  f  after  running  the  Levenberg-Marquardt  routine  using  a  frequency  of  10 6 Hz.  Note  that  while  the 
initial  values  in  case  1  are  farthest  from  the  exact  solution,  some  of  the  corresponding  t  are  actually  closest 
to  the  actual  value  (for  example,  oq  =  .1).  While  this  suggests  that  this  should  more  easily  converge  to  the 
“curve  of  best  fit” ,  it  will  most  likely  be  farther  away  from  the  exact  solution  on  this  curve,  and  therefore 
should  still  be  considered  the  hardest  of  the  cases  to  solve. 

Recognizing  that  we  have  converged  to  the  “curve  of  best  fit”  in  each  of  the  above  cases,  we  may  now 
restart  an  optimization  routine  that  traverses  this  curve.  Some  modifications  to  the  optimization  routine’s 
parameters  are  required  to  address  the  vast  difference  in  scales  of  this  subproblem  as  compared  to  the  two 
parameter  optimization  problem.  The  final  results  of  this  two  step  optimization  approach  are  given  in  Tables 
13  through  18.  As  expected,  the  results  from  case  3  are  generally  better  than  the  other  two  sets  of  initial 
conditions.  The  highest  frequency  attempted,  \QllHz  seemed  to  perform  the  most  poorly.  This  suggests 
that  the  higher  the  frequency,  the  more  difficult  to  accurately  resolve  the  polarization  mechanisms.  Although 
the  106Hz  case  used  the  curve  of  constant  t  while  the  10 9Hz  case  used  the  curve  of  constant  A,  there  was 
no  evidence  to  suggest  that  one  case  performed  better  than  the  other.  Lastly,  it  appears  that  when  material 
1  (corresponding  to  relaxation  parameter  T\)  is  of  the  highest  proportion,  the  optimization  routine  is  best 
able  to  resolve  iq.  Likewise,  if  material  1  is  of  a  lower  proportion,  the  routine  instead  does  a  better  job  of 
resolving  t^.  Note  that  there  are  several  instances  where  the  optimizer  “switched”  iq  and  72,  for  example  in 
case  1  of  frequency  10 9 Hz  with  oq  =  .5.  In  this  scenario  each  material  comprises  50%  of  the  whole,  so  the 
problem  is  symmetric  and  swapping  n  and  Tq  has  no  effect.  However,  in  case  1  of  frequency  10 9 Hz  with 
oq  =  .1,  it  appears  that  t\  is  converging  toward  the  exact  t%  value,  but  since  this  problem  is  not  symmetric, 
72  converges  to  a  meaningless  solution. 
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Table  8:  The  initial  values  of  A  :=  ^  +  ^jr  for  each  set  of  initial  conditions  (case  0  represents  the  exact 
solution) . 


case 

.1 

.5 

.9 

0 

0.200556 

0.158333 

0.116111 

1 

0.952112 

0.538334 

0.124556 

2 

0.385278 

0.237500 

0.0897223 

3 

0.245945 

0.174167 

0.102389 

Table  9:  Resulting  values  of  A  from  the  Levenberg-Marquardt  routine  using  a  frequency  of  10 11  Hz  for  each 
set  of  initial  conditions  (case  0  represents  the  exact  solution) . 


case 

.1 

d  i 

.5 

.9 

0 

0.200556 

0.158333 

0.116111 

1 

0.200573 

0.15834 

0.116109 

2 

0.200573 

0.158327 

0.1159 

3 

0.200573 

0.158363 

0.116114 

Table  10:  Resulting  values  of  A  from  the  Levenberg-Marquardt  routine  using  a  frequency  of  10 9 Hz  for  each 
set  of  initial  conditions  (case  0  represents  the  exact  solution) . 


case 

.1 

d  i 

.5 

.9 

0 

0.200556 

0.158333 

0.116111 

1 

0.200555 

0.158331 

0.116029 

2 

0.200556 

0.158337 

0.116132 

3 

0.200556 

0.158339 

0.116119 

Table  11:  The  initial  values  of  f  :=  ol\T\  +  a^Tn  for  each  set  of  initial  conditions  (case  0  represents  the  exact 
solution) . 


case 

.1 

.5 

.9 

0 

1.7380e-08 

2.3700e-08 

3.0020e-08 

1 

1.8644e-08 

8.0580e-08 

1.42516e-07 

2 

1.3430e-08 

3.5550e-08 

5.7670e-08 

3 

1.5326e-08 

2.6070e-08 

3.6814e-08 

Table  12:  Resulting  values  of  f  from  the  Levenberg-Marquardt  routine  using  a  frequency  of  10 6Hz  for  each 
set  of  initial  conditions  (case  0  represents  the  exact  solution) . 


case 

.1 

d\ 

.5 

.9 

0 

1.7380e-08 

2.3700e-08 

3.0020e-08 

1 

1.73954e-08 

2.37014e-08 

3.00956e-08 

2 

1.73819e-08 

2.37031e-08 

3.00222e-08 

3 

1.73803e-08 

2.37016e-08 

3.00241e-08 
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Table  13:  Final  estimates  for  t\  from  two  step  optimization  approach  using  a  frequency  of  10 llHz  for  each 
set  of  initial  conditions  (recall  the  exact  solution  T-f  =3.1600e-8). 


case 

.1 

d 1 

.5 

.9 

i 

1.57221e-07 

2.81579e-08 

2.61108e-08 

2 

6.34951e-08 

4.63981e-08 

3.15232e-08 

3 

3.98043e-08 

3.33803e-08 

3.2833e-08 

Table  14:  Final  estimates  for  n  from  two  step  optimization  approach  using  a  frequency  of  109Hz  for  each 
set  of  initial  conditions  (recall  the  exact  solution  =3.1600e-8). 


case 

.1 

d\ 

.5 

.9 

i 

1.12709e-08 

1.58005e-08 

3.15929e-08 

2 

3.1574e-08 

3.16002e-08 

3.16017e-08 

3 

3.16009e-08 

3.16009e-08 

3.16009e-08 

Table  15:  Final  estimates  for  t\  from  two  step  optimization  approach  using  a  frequency  of  10 eHz  for  each 
set  of  initial  conditions  (recall  the  exact  solution  T-f  =3.1600e-8). 


case 

.1 

d\ 

.5 

.9 

i 

3.17673e-08 

3.16031e-08 

3.16986e-08 

2 

3.16206e-08 

3.16068e-08 

3.16031e-08 

3 

3.16031e-08 

3.16039e-08 

3.16053e-08 

Table  16:  Final  estimates  for  r2  from  two  step  optimization  approach  using  a  frequency  of  10 llHz  for  each 
set  of  initial  conditions  (recall  the  exact  solution  r|  =1.5800e-8). 


case 

.1 

d  i 

.5 

.9 

i 

1.51273e-08 

1.68275e-08 

2.93245e-07 

2 

1.53699e-08 

1.36276e-08 

1.61391e-08 

3 

1.56196e-08 

1.53854e-08 

1.35142e-08 

Table  17:  Final  estimates  for  r2  from  two  step  optimization  approach  using  a  frequency  of  10 9 Hz  for  each 
set  of  initial  conditions  (recall  the  exact  solution  r|  =1.5800e-8). 


case 

.1 

a  i 

.5 

.9 

i 

1.75594e-08 

3.15995e-08 

1.58782e-08 

2 

1.58008e-08 

1.57994e-08 

1.57801e-08 

3 

1.58001e-08 

1.5799e-08 

1.57925e-08 

Table  18:  Final  estimates  for  r2  from  two  step  optimization  approach  using  a  frequency  of  10 6Hz  for  each 
set  of  initial  conditions  (recall  the  exact  solution  r|  =1.5800e-8). 


case 

.1 

d  i 

.5 

.9 

i 

1.52384e-08 

1.52286e-08 

1.61666e-08 

2 

1.5787e-08 

1.45007e-08 

1.18861e-08 

3 

1.57845e-08 

1.55744e-08 

1.5783e-08 
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5.3  Determination  of  Volume  Proportions 

We  now  attempt  to  apply  a  one  parameter  Levenberg-Marquardt  optimization  routine  to  our  problem  to 
identify  the  relative  amounts  of  two  materials  with  known,  distinct  relaxation  times.  Thus  we  are  trying  to 
find  the  corresponding  volume  proportions  of  the  two  materials  and  ai  :=  1  —  a\).  We  again  consider  the 
three  scenarios  with  respect  to  the  exact  volume  proportions:  ot\  £  {.1,  .5,  .9}.  We  also  perform  our  inverse 
problem  using  the  frequencies  1011  Hz,  10 9 Hz  and  10 6 Hz.  Lastly  we  test  the  optimization  procedure  with 
three  various  initial  conditions:  a3  G  {.9999,  .0001,  .5}  (except  in  the  case  when  a\  =  .5  in  which  case  we 
used  a®  £  {.9999,  .0001,  .4}).  We  refer  to  these  as  Cases  1,  2,  and  3,  respectively.  In  all  of  the  following  we 
assume  that  the  known  relaxation  times  are  t\  =  io_7-50031  ss  3.16  x  10-8  and  T2  =  io-7-80134  ss  1.58  x  10-8. 

Figures  11  through  13  depict  the  graphs  of  the  functions  that  we  are  attempting  to  minimize.  For  each 
case  the  curves  appear  well  behaved.  The  results  for  this  one  parameter  inverse  problem,  displayed  in  Tables 
19  through  21,  verify  that  the  relative  proportions  of  known  materials  are  generally  easily  identifiable.  Tables 
22  through  24  display  the  final  objective  function  values  for  each  case.  Note  that  typical  initial  values  for  J 
were  around  0.1,  further,  the  tolerance  was  set  at  10~9,  thus  all  but  a  few  cases  converged  before  reaching 
the  maximum  of  20  iterations. 


f=1e11,a1=.1 


Figure  11:  The  objective  function  for  the  relaxation  time  inverse  problem  versus  oq  using  a  frequency  of 
\QllHz  and  =  .1. 
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0.025 


f=1e11,a  =.5 


Figure  12:  The  objective  function  for  the  relaxation  time  inverse  problem  versus  oq  using  a  frequency  of 
10 lxHz  and  aj  =  .5. 


f=1e1 1  ,ct  =.9 


Figure  13:  The  objective  function  for  the  relaxation  time  inverse  problem  versus  oq  using  a  frequency  of 
\QllHz  and  =  .9. 
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Table  19:  Results  for  the  one  parameter  inverse  problem  to  determine  the  relative  proportion  of  two  known 
Debye  materials  using  a  frequency  of  10 11  Hz  (aq  estimates  are  shown). 


case 

.1 

d\ 

.5 

.9 

i 

0.10095 

0.501353 

0.900013 

2 

0.0999994 

0.5 

0.9 

3 

0.100643 

0.499994 

0.899994 

Table  20:  Results  for  the  one  parameter  inverse  problem  to  determine  the  relative  proportion  of  two  known 
Debye  materials  using  a  frequency  of  10 9Hz  (aq  estimates  are  shown). 


case 

.1 

d\ 

.5 

.9 

i 

0.10086 

0.5 

0.9 

2 

0.0999161 

0.5 

0.9 

3 

0.100017 

0.5 

0.9 

Table  21:  Results  for  the  one  parameter  inverse  problem  to  determine  the  relative  proportion  of  two  known 
Debye  materials  using  a  frequency  of  10 6Hz  ( oq  estimates  are  shown). 


case 

.1 

d  i 

.5 

.9 

i 

0.0995386 

0.5 

0.900008 

2 

0.0999965 

0.499969 

0.899996 

3 

0.100013 

0.499999 

0.899999 

Table  22:  Final  objective  function  values  for  the  inverse  problem  to  determine  the  relative  proportion  of  two 
known  Debye  materials  using  a  frequency  of  \QllHz. 


case 

.1 

d  i 

.5 

.9 

i 

8.27897e-08 

1.69802e-07 

1.81577e-ll 

2 

3.46217e-14 

1.64325e-18 

1.8982e-15 

3 

3.79016e-08 

3.76649e-12 

2.91698e-12 

Table  23:  Final  objective  function  values  for  the  inverse  problem  to  determine  the  relative  proportion  of  two 
known  Debye  materials  using  a  frequency  of  10 9 Hz. 


case 

.1 

d\ 

.5 

.9 

i 

1.63248e-05 

1.25127e-12 

6.89036e-14 

2 

1.55308e-07 

5.89446e-13 

5.8576e-14 

3 

6.53272e-09 

2.37274e-12 

5.77714e-15 

Table  24:  Final  objective  function  values  for  the  inverse  problem  to  determine  the  relative  proportion  of  two 
known  Debye  materials  using  a  frequency  of  10 eHz. 


case 

.1 

d\ 

.5 

.9 

i 

7.74602e-08 

7.31725e-15 

4.39638e-12 

2 

4.51572e-12 

1.40031e-10 

1.03552e-12 

3 

6.11308e-ll 

2.0608e-13 

1.53426e-13 
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5.4  Determination  of  Volume  Proportions  and  Relaxation  Times  Simultane¬ 
ously 

Although  we  anticipate  this  problem  formulation  to  be  underdetermined,  we  run  our  optimization  routine  for 
the  problem  where  neither  the  relaxation  times,  nor  the  relative  volume  proportions,  of  two  distinct  Debye 
materials  are  known.  Thus  this  is  a  three  parameter  inverse  problem  for  {ti,t2,o:i}  (since  a2  =  1  —  aq). 
Again  we  consider  the  following  scenarios  with  respect  to  the  actual  volume  proportions:  aq  £  {.1,  .5,  .9}. 
We  also  perform  our  inverse  problem  using  the  frequencies  10 llHz,  109Hz  and  10 6Hz.  Lastly  we  test  the 
optimization  procedure  with  the  three  various  initial  conditions  that  were  given  in  Table  1.  Note  that  the 
actual  values  are  still  n  =  lO"7-50031  «  3.16  x  10"8  and  r2  =  lO"7-80134  «  1.58  x  10"8. 

Our  initial  condition  for  the  volume  distribution  is  a3  =  .9999,  which  means  there  is  essentially  only 
material  1.  We  chose  this  initial  condition  mainly  to  see  whether  the  optimization  routine  would  try  to  com¬ 
pensate  for  the  altered  volume  distribution  by  significantly  changing  the  relaxation  times  and  not  sufficiently 
correcting  the  volume  distribution.  This  is  in  fact  what  seems  to  have  occurred.  The  final  aq  values  are 
given  in  Tables  25  through  27,  while  the  final  T\  and  t2  values  are  in  Tables  28  through  30  and  31  through 
33  respectively.  Tables  34  through  36  display  the  final  objective  function  values.  Any  values  of  the  final 
objective  function  under  the  tolerance  of  10-9  should  be  considered  as  indicative  of  convergence,  while  values 
greater  than  1  clearly  indicate  stagnation,  (e.g.,  10 9Hz:  case  1). 


Table  25:  Resulting  values  of  aq  for  the  underdetermined  inverse  problem  using  a  frequency  of  lQllHz. 


case 

.1 

a\ 

.5 

.9 

i 

0.606255 

0.881722 

0.912728 

2 

0.761665 

0.817744 

0.877667 

3 

0.818091 

0.879896 

0.942175 

Table  26:  Resulting  values  of  oq  for  the  underdetermined  inverse  problem  using  a  frequency  of  10 9 Hz. 


case 

.1 

.5 

.9 

i 

0.941041 

0.941091 

0.880293 

2 

0.999811 

0.993342 

0.994368 

3 

0.990988 

0.96514 

0.965754 

Table  27:  Resulting  values  of  oq  for  the  underdetermined  inverse  problem  using  a  frequency  of  10 6Hz. 


case 

.1 

CH  i 

.5 

.9 

i 

0.687959 

0.746478 

0.956578 

2 

0.83678 

0.900072 

0.908588 

3 

0.890569 

0.841634 

0.945684 

Of  the  cases  that  converged,  none  converged  to  the  correct  solution.  For  example,  the  final  values  of  aq 
were  all  greater  than  .6,  (although  the  estimates  corresponding  to  a  smaller  aj  were  on  average  less  than  the 
estimates  corresponding  to  aj  =  .9).  When  looking  at  the  final  estimates  for  the  relaxation  times,  however, 
there  appears  to  be  no  rhyme  nor  reason.  This  is  quite  similar  to  what  happened  when  we  first  tried  to 
determine  unknown  relaxation  times,  but  with  a  fixed  volume  proportion,  in  Section  5.  Thus  we  may  expect 
that  here  again  no  single  parameter,  n,  r2,  nor  aq,  may  converge,  but  instead  we  may  see  convergence  of 
a  general  relation  involving  them  all,  namely  A  :=  ^  (or  respectively  f  :=  a,\T\  +  a2r2  for  angular 

frequencies  less  than  4). 

We  redisplay  the  exact  values  of  A  and  f  in  Table  37  for  reference.  Also,  the  initial  values  of  A  and  f 
for  each  case  of  initial  conditions  (recall  Table  1)  is  given  in  Table  38.  Finally,  Tables  39  through  41  display 
the  error  of  the  final  resulting  A  (or  f  appropriately)  estimates  from  the  exact  solutions  (i.e. ,  the  absolute 
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Table  28:  Resulting  values  of  iq  for  the  underdetermined  inverse  problem  using  a  frequency  of  \QllHz  (recall 
the  exact  solution  r*  =3.1600e-8). 


case 

.1 

aii 

.5 

.9 

i 

3.9469e-08 

5.28765e-07 

3.23635e-08 

2 

1.52948e-08 

2.49314e-08 

2.60749e-08 

3 

1.79134e-08 

2.40213e-08 

3.06153e-08 

Table  29:  Resulting  values  of  T\  for  the  underdetermined  inverse  problem  using  a  frequency  of  10 9 Hz  (recall 
the  exact  solution  r*  =3.1600e-8). 


case 

.1 

ai 

.5 

.9 

i 

1.2046e-07 

1.24415e-07 

5.50848e-07 

2 

1.67567e-08 

2.16631e-08 

2.93323e-08 

3 

1.65557e-08 

2.25001e-08 

3.0285e-08 

Table  30:  Resulting  values  of  iq  for  the  underdetermined  inverse  problem  using  a  frequency  of  10 6Hz  (recall 
the  exact  solution  t{  =3.1600e-8). 


case 

.1 

ai 

.5 

.9 

i 

2.12522e-08 

2.81013e-08 

3.10195e-08 

2 

1.91221e-08 

2.3216e-08 

3.15199e-08 

3 

2.58248e-08 

2.67889e-08 

3.11248e-08 

values  of  the  differences).  We  show  the  errors  instead  of  the  actual  values  because  some  cases  converged 
so  well  that  the  number  of  decimal  places  required  to  see  discrepancies  is  impractical  to  show!  Thus  while 
the  inverse  problem  was  unable  to  accurately  resolve  the  individual  values  iq,  T2,  and  a\,  in  all  but  a  few 
cases  (particularly  in  the  10 9 Hz  scenario)  the  optimization  routine  did  converge  as  well  as  to  be  expected  to 
the  “curve  of  best  fit”.  In  the  previous  section,  with  the  fixed  volume  proportions,  we  were  able  to  traverse 
the  curve  of  constant  A  (or  t  appropriately)  to  find  the  global  minimum.  However,  here  we  truly  have  an 
underdetermined  problem  and  therefore  we  cannot  extract  any  further  information  without  providing  more 
data. 
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Table  31:  Resulting  values  of  72  for  the  underdetermined  inverse  problem  using  a  frequency  of  10 11  Hz  (recall 
the  exact  solution  r|  =1.5800e-8). 


case 

.1 

ai 

.5 

.9 

1 

8.79556e-09 

2.58245e-09 

1.32076e-08 

2 

2.30775e-08 

1.24249e-08 

1.06318e-07 

3 

1.25825e-08 

1.10813e-08 

1.4329e-08 

Table  32:  Resulting  values  of  72  for  the  underdetermined  inverse  problem  using  a  frequency  of  10 9 Hz  (recall 
the  exact  solution  =1.5800e-8). 


case 

.1 

a  1 

.5 

.9 

1 

1.81687e-10 

2.32623e-10 

1.2353e-09 

2 

4.08525e-10 

4.11677e-09 

6.18566e-09 

3 

3.36629e-08 

7.61993e-09 

1.17231e-08 

Table  33:  Resulting  values  of  72  for  the  underdetermined  inverse  problem  using  a  frequency  of  10 6Hz  (recall 
the  exact  solution  r|  =1.5800e-8). 


case 

.1 

.5 

.9 

1 

8.82441e-09 

1.06962e-08 

7.87058e-09 

2 

8.40313e-09 

2.78555e-08 

1.5107e-08 

3 

6.19113e-09 

7.15491e-09 

1.0713e-08 

Table  34:  Resulting  values  of  the  objective  function  J  for  the  underdetermined  inverse  problem  using  a 
frequency  of  10 xlHz. 


case 

.1 

a\ 

.5 

.9 

1 

1.87304e-09 

1.39638e-07 

2.96012e-12 

2 

1.87016e-15 

3.3089e-14 

1.60601e-17 

3 

2.79754e-14 

2.80156e-16 

1.38435e-13 

Table  35:  Resulting  values  of  the  objective  function  J  for  the  underdetermined  inverse  problem  using  a 
frequency  of  10 9 Hz. 


case 

.1 

.5 

.9 

1 

9.59505 

13.787 

25.045 

2 

0.00021558 

2.49585e-06 

2.03659e-07 

3 

3.48546e-06 

5.0347e-07 

1.18429e-08 

Table  36:  Resulting  values  of  the  objective  function  J  for  the  underdetermined  inverse  problem  using  a 
frequency  of  10 6 Hz. 


case 

.1 

.5 

.9 

1 

1.19225e-05 

4.76542e-07 

1.18516e-07 

2 

6.34493e-07 

3.86474e-05 

4.41757e-10 

3 

4.07294e-06 

1.99923e-06 

2.75692e-08 
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Table  37:  The  exact  values  of  A  :=  ^  ^  (first  row)  and  r  :=  aq n  +  (second  row)  for  each  set  of 

volume  distributions.  _ 


* 

a  i 

.1 

.5 

.9 

A* 

0.200556 

0.158333 

0.116111 

1.7380e-08 

2.3700e-08 

3.0020e-08 

Table  38:  The  initial  values  of  A  :=  ^  ^  (first  column)  and  t  :=  ol\T\  +  a2T2  (second  column)  for  each 

set  of  initial  conditions.  _ 


case 

A0 

f° 

i 

0.0212146 

1.57985e-07 

2 

0.0528147 

6.31945e-08 

3 

0.0844624 

3.94973e-08 

Table  39:  Error  of  the  resulting  values  of  A  from  the  exact  values  for  the  underdetermined  inverse  problem 
using  a  frequency  of  10 11  Hz  for  each  set  of  initial  conditions. 


case 

.1 

a  i 

.5 

.9 

i 

2.99851e-08 

7.35704e-08 

3.49983e-09 

2 

1.23759e-10 

4.68118e-10 

2.30235e-10 

3 

1.14292e-09 

2.16657e-10 

1.273e-10 

Table  40:  Error  of  the  resulting  values  of  A  from  the  exact  values  for  the  under  determined  inverse  problem 

lioirwr  o  fronnonmr  rvf  1  Ur/  fnr  oonVi  oot  rq-p  initial  pati  pi  i  T  i  An  c: 


case 

.1 

o/.  i 

.5 

.9 

i 

0.907928 

0.711589 

0.212452 

2 

4.93369e-06 

1.05337e-05 

1.47737e-06 

3 

3.86466e-06 

4.68234e-06 

6.40273e-08 

Table  41:  Error  of  the  resulting  values  of  f  from  the  exact  values  for  the  underdetermined  inverse  problem 

lioirwr  o  frnrmnnmr  rvf  1  fnr  oonVi  oof  rvf  initial  PAn pi  i  T  i An o 


case 

.1 

d\ 

.5 

.9 

i 

5.78065e-12 

1.12783e-ll 

5.70111e-12 

2 

7.42983e-12 

2.03867e-ll 

4.54227e-13 

3 

2.3737e-ll 

2.04292e-ll 

3.88617e-12 
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5.5  Multiple  Frequency  Interrogation 

In  certain  situations  it  may  be  feasible  to  add  to  the  data  used  in  the  inverse  problem  in  the  previous 
section  by  interrogating  the  material  again  with  a  different  frequency.  In  particular,  if  for  example  the  first 
interrogating  angular  frequency  is  less  than  c oc  :=  1/r  then  choosing  the  second  frequency  to  be  greater  than 
this  critical  value  should  give  different  information,  i.e. ,  lol  <  1/r  <  ojh-  We  have  already  seen  that  for 
frequencies  lower  than  the  critical  angular  frequency  we  are  able  to  accurately  resolve  f.  Thus,  using  this 
information  we  may  choose  a  second  frequency  to  be  higher  than  the  critical  level,  and  therefore  expect  to 
be  able  to  resolve  A,  as  was  done  in  our  previous  examples  for  high  frequencies.  We  applied  this  approach 
to  the  test  cases  described  for  the  /  =  106 Hz  interrogating  signal.  Using  the  results  of  the  previous  section, 
i.e.,  fixing  f  and  using  the  corresponding  values  of  ti,  72,  and  oq  as  initial  conditions,  we  attempted  an 
inverse  problem  to  match  data  obtained  by  interrogating  the  same  material  with  a  frequency  /  =  109  Hz, 
which  is  higher  than  the  critical  value.  In  particular,  we  performed  a  one  parameter  search  for  each  of  T\  and 
r 2  successively,  then  a  two  parameter  search  for  both  simultaneously  for  fine-tuning.  As  in  previous  inverse 
problems,  this  allowed  determination  of  A.  Now  with  A  and  t  both  fixed,  we  traverse  along  both  contours 
simultaneously  with  a  one  parameter  inverse  problem  for  t\.  Lastly,  again,  we  do  a  full  three  parameter 
inverse  problem  for  fine-tuning.  The  results  for  Cases  1,  2,  and  3  are  given  in  Tables  42  -  44.  Note  that  in 
order  to  ensure  better  accuracy  when  the  wavelength  is  decreased,  we  double  the  number  of  finite  elements 
used  to  solve  the  high  frequency  inverse  problem. 

In  each  case  the  lower  frequency  inverse  problem  was  able  to  determine  the  value  of  f  to  at  least  three 
decimal  places,  while  the  values  of  t\,  72,  oq  and  A  are  in  general  not  necessarily  even  improved  over  the 
initial  estimates  (e.g.,  T\  for  Case  3).  This  is  similar  to  the  results  of  previous  sections  before  traversing  the 
contours  of  constant  f  or  A.  The  higher  frequency  inverse  problem  with  t  fixed  was  able  to  determine  the 
value  of  A  to  at  least  three  decimal  places  in  each  case,  while  again  the  values  of  T\,  T2  and  aq  are  in  general 
not  necessarily  even  improved  over  the  previous  estimates.  The  final  rows  show  the  results  from  running  a 
one  parameter  inverse  problem  for  t\  with  both  f  and  A  fixed  (and  then  using  a  three  parameter  search  for 
fine  tuning).  Each  of  the  final  estimates  for  rq,  72,  and  oq  have  at  least  two  decimal  places  accuracy. 

Recall  that,  as  explained  in  Remark  1,  the  slab  thickness  here  is  400m.  These  dimensions  are  merely 
arbitrary  examples,  but  there  is  a  direct  restriction  on  the  width  of  a  slab  if  it  is  to  be  reliably  simultaneously 
interrogated  with  frequencies  just  above  and  just  below  the  critical  angular  frequency  for  a  given  medium, 
which  is  described  by  the  material  properties  through  the  value  of  f.  Thus,  for  example,  materials  with 
relaxation  times  on  the  order  of  10-11,  like  water,  that  have  dimensions  on  the  order  of  .4m  are  quite 
feasible  to  interrogate  with  multiple  frequencies  above  and  below  its  critical  frequency. 
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Table  42:  Numerical  results  from  the  inverse  problem  using  data  from  interrogating  the  same  medium  with 
one  frequency  above,  and  one  frequency  below  the  critical  angular  frequency  for  that  medium.  Given  are 
the  ti,  r 2,  and  aq  values  along  with  the  relative  error  for  A  and  f  for  each  step  of  the  inverse  problem:  exact 
solution,  initial  conditions  for  Case  1,  and  the  solutions  from  the  low  frequency  inverse  problem,  from  the 
high  frequency  inverse  problem  with  fixed  f ,  and  finally  from  fixing  both  f  and  A. 


Case  1 

T\ 

T2 

Oil 

relerr(A) 

relerr(f ) 

Exact 

3.16e-8 

1.58e-8 

0.1 

0 

0 

Initial 

1.58e-7 

3.16e-9 

0.9999 

0.894221 

8.09005 

/  =  10  6Hz 

2.12522e-8 

8.82441e-9 

0.687957 

0.126493 

3.33717e-4 

f  =  10 9  Hz 

1.82157e-8 

8.82198e-9 

0.910418 

1.24653e-4 

3.33717e-4 

Contour 

3.13819e-8 

1.57826e-8 

0.102596 

4.98614e-5 

1.72612e-4 

Table  43:  Numerical  results  from  the  inverse  problem  using  data  from  interrogating  the  same  medium  with 
one  frequency  above,  and  one  frequency  below  the  critical  angular  frequency  for  that  medium.  Given  are 
the  ti,  T2,  and  op  values  along  with  the  relative  error  for  A  and  f  for  each  step  of  the  inverse  problem:  exact 
solution,  initial  conditions  for  Case  2,  and  the  solutions  from  the  low  frequency  inverse  problem,  from  the 
high  frequency  inverse  problem  with  fixed  r,  and  finally  from  fixing  both  t  and  A. 


Case  2 

T\ 

t2 

Oil 

relerr(A) 

relerr(f) 

Exact 

3.16e-8 

1.58e-8 

0.1 

0 

0 

Initial 

6.32e-8 

7.9e-9 

0.9999 

0.736658 

2.63605 

/  =  10 6  Hz 

1.91221e-8 

8.40313e-9 

0.836785 

5.08337e-2 

4.25777e-4 

f  =  10 9  Hz 

1.81458e-8 

8.47912e-9 

0.920012 

1.2964e-4 

4.25777e-4 

Contour 

3.11405e-8 

1.57683e-8 

0.105044 

5.48475e-5 

1.78366e-4 

Table  44:  Numerical  results  from  the  inverse  problem  using  data  from  interrogating  the  same  medium  with 
one  frequency  above,  and  one  frequency  below  the  critical  angular  frequency  for  that  medium.  Given  are 
the  ti,  T2,  and  oq  values  along  with  the  relative  error  for  A  and  t  for  each  step  of  the  inverse  problem:  exact 
solution,  initial  conditions  for  Case  3,  and  the  solutions  from  the  low  frequency  inverse  problem,  from  the 
high  frequency  inverse  problem  with  fixed  r,  and  finally  from  fixing  both  f  and  A. 


Case  3 

Tl 

T2 

Oil 

relerr(A) 

relerr(f ) 

Exact 

3.16e-8 

1.58e-8 

0.1 

0 

0 

Initial 

3.95e-8 

1.264e-9 

0.9999 

0.578859 

1.27257 

/  =  10 6  Hz 

1.81785e-8 

1.09823e-8 

0.887927 

1.79102e-2 

4.60299e-4 

f  =  10 9  Hz 

1.87871e-8 

1.09078e-8 

0.820402 

1.09695e-4 

4.60299e-4 

Contour 

3.12584e-8 

1.57761e-8 

0.103767 

4.48752e-5 

1.49597e-4 
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6  Inverse  Problem  Results  Using  a  Uniform  Distribution 

The  previous  polarization  inverse  problems  have  assumed  a  discrete  distribution  with  two  atoms.  Accord¬ 
ing  to  experimental  reports  [23],  most  materials  demonstrate  polarization  effects  described  by  a  range  of 
relaxation  times.  Here  we  consider  the  simplest  of  distributions  by  exploring  the  possibility  of  a  uniform 
distribution  of  relaxations  times  (r)  between  a  lower  and  upper  limit  (ra  and  77,  respectively).  This  presents 
us  again  with  a  two  parameter  inverse  problem,  namely,  to  try  to  resolve  the  end  points  of  the  distribution 
used  to  generate  the  given  data.  Computationally  we  still  approximate  this  distribution  with  discrete  nodes, 
but  instead  of  just  one  at  each  endpoint,  we  use  l  =  13  uniformly  distributed  within  the  interval  (recall 
Section  3.2).  Note  that  we  do  not  wish  to  restrict  our  optimization  routine  to  search  only  for  ra  <  17,, 
therefore  subsequently  if  ra  >  77,  it  is  assumed  without  loss  of  generality  that  77,  is  the  lower  limit  of  the 
distribution  and  ra  is  the  upper  limit. 

Figures  14  and  15  depict  the  objective  function  and  the  log  of  the  objective  function,  respectively,  for  a 
frequency  of  1( l11  Hz.  The  solid  line  in  Figure  15  is  the  curve  of  constant  A.  Given  a  uniform  distribution  of 
relaxation  times  in  a  Debye  medium,  this  parameter  is  given  (analytically)  by 

r  -=  1  fTb  dr  =  In  77,  -  In  ra 

c{n  -  To)  JTa  T  c(rb  -  To) 

Note  that  in  computations  we  must  use  the  same  quadrature  method  to  evalute  this  integral  as  we  do  for  Q 
to  ensure  a  correct  correlation.  Although  the  curve  appears  slightly  different  from  the  discrete  distribution 
case  in  Figure  3,  the  fact  that  this  objective  function  is  also  small  along  this  curve  suggests  that  this  problem 
should  behave  similarly  to  the  discrete  distribution  case  in  Section  5.  Figures  16  and  17  depict  the  objective 
function  and  the  log  of  the  objective  function,  respectively,  for  a  frequency  of  10 6Hz.  We  notice  that  again 
the  orientation  of  the  “curve  of  best”  fit  is  different  from  the  higher  frequency  case.  The  solid  line  in  Figure 
17  is  the  curve  of  constant  f  :=  f^b  rdF(r )  =  Tbl2Ta  ■  Again,  the  fact  that  the  A  term  only  dominates  the 
behavior  when  the  interrogating  frequency  is  greater  than  is  consistent  with  the  discrete  distribution 
case. 

Based  on  our  previous  experience  with  the  discrete  distribution,  we  anticipate  that  the  two  parameter 
inverse  problem  will  simply  converge  to  the  “line  of  best  fit”.  Thus  instead,  we  use  an  equivalent  method  to 
converge  to  this  curve,  namely  minimizing  over  just  one  parameter,  ra,  leaving  rb  fixed  at  its  initial  value  (the 
same  values  as  those  given  as  r9  in  Table  1).  It  should  be  noted  that  one  may  actually  perform  two  of  the 
one  parameter  constrained  optimizations,  one  in  each  of  the  directions  ra  and  rb,  to  allow  for  the  possibility 
of  the  first  direction  not  converging,  for  example,  if  t9,t9  <  ICC9.  In  general,  the  second  direction  of  the 
optimization  requires  only  enough  iterations  to  verify  convergence,  as  the  first  direction  has  usually  already 
converged  to  the  “line  of  best  fit”.  Using  this  one  parameter  at  a  time  approach,  we  still  converge  to  the  “line 
of  best  fit” ,  but  theoretically  use  half  as  many  function  evaluations  as  the  two  parameter  inverse  problem, 
since  we  only  compute  one  gradient  at  each  step.  (In  practice,  only  a  third  as  many  function  evaluations 
were  actually  needed  to  get  the  same  order  of  accuracy  as  the  two  parameter  inverse  problem.) 

We  performed  the  one  parameter  inverse  problem  using  the  three  initial  condition  cases  described  above 
in  Table  1  and  again  using  frequencies  lO^Hz,  10 ^ Hz  and  10 °Hz.  The  r0  estimates  from  running  Levenberg- 
Marquardt  on  the  modified  least  squares  objective  function  are  given  in  Table  45.  Again,  as  in  the  discrete 
case,  the  values  of  the  relaxation  times  do  not  appear  to  be  converging  to  the  correct  solution.  But  we 
expect  that  the  optimization  routine  is  converging  to  the  “curve  of  best  fit” .  To  test  this  we  must  look  at 
the  approximations  to  A  and  f . 

The  initial  values  of  A  and  t  are  given  in  Table  46  (note  that  these  values  were  computed  by  A  «  ]C,-  — 
and  t  «  cnn  using  appropriately  defined  determined  by  the  Composite  Simpson’s  rule).  The 

exact  values  of  each  are  A*  =  0.248369  and  f*  =  5.42467  x  10-8.  The  resulting  A  and  f  values  from  running 
the  one  parameter  Levenberg-Marquardt  routine  are  given  in  Table  47.  Clearly  each  case  has  converged  to 
the  “line  of  best  fit” ;  in  general  the  closer  initial  conditions  converged  closer  to  the  actual  value  of  A  (or  t 
for  /  =  10 6 Hz). 

After  our  one  parameter  optimization  routine  resolved  A  (or  t  for  /  =  106  Hz),  we  minimized  for  r0  along 
the  line  of  constant  A  (or  f).  Again,  this  is  a  one  parameter  inverse  problem,  and  therefore  very  efficient. 
The  results  of  these  computations  are  given  in  Tables  48  and  49,  for  r0  and  the  corresponding  rb  (given 
constant  A  or  f),  respectively. 
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Figure  15:  The  log  of  the  objective  function  for  the  uniform  distribution  inverse  problem  versus  the  log  of 
Ta  and  the  log  of  t>,  using  a  frequency  of  10nHz.  The  solid  line  above  the  surface  represents  the  curve  of 
constant  A. 


iog(taub) 

Figure  14:  The  objective  function  for  the  uniform  distribution  inverse  problem  versus  the  log  of  ra  and  the 
log  of  Tb  using  a  frequency  of  10nHz. 
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J  with  block  distribution  (10A-7.5,10A-7.8)  using  f=1e6 
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Figure  16:  The  objective  function  for  the  uniform  distribution  inverse  problem  versus  the  log  of  ra  and  the 
log  of  Tb  using  a  frequency  of  106iJ^. 


log(J)  with  block  distribution  (10A-7.5,10A-7.8)  using  f=1e6 


Figure  17:  The  log  of  the  objective  function  for  the  uniform  distribution  inverse  problem  versus  the  log  of 
Ta  and  the  log  of  Tb  using  a  frequency  of  106 Hz.  The  solid  line  above  the  surface  represents  the  curve  of 
constant  t . 
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Table  45:  Resulting  values  of  rQ  from  the  one  parameter  Levenberg-Marquardt  routine  for  the  inverse 
problem  to  determine  the  endpoints  of  a  uniform  distribution  of  relaxation  times  (recall  the  exact  solution 
r*  =3.16000e-8).  _ 


case 

Initial 

1011 

Frequency  (Hz) 
109 

106 

i 

1.58000e-7 

5.41874e-8 

5.40910e-8 

3.59781e-8 

2 

6.32000e-8 

4.01846e-8 

4.01788e-8 

3.43510e-8 

3 

3.95000e-8 

3.41828e-8 

3.41820e-8 

3.27009e-8 

Table  46:  The  initial  values  of  A  :=  JT  —  (or  f  :=  JT  aiTi  for  /  =  106)  for  each  set  of  initial  conditions 
for  the  inverse  problem  to  determine  the  endpoints  of  a  uniform  distribution  of  relaxation  times  (case  0 
represents  the  exact  solution). 


case 

Frequency  (Hz) 
1011  -  109 

106 

0 

0.248369 

5.42467e-8 

1 

0.115739 

2.33313e-7 

2 

0.176705 

9.66433e-8 

3 

0.223136 

6.42533e-8 

Table  47:  Resulting  values  of  A  (or  f  for  /  =  106)  from  the  Levenberg-Marquardt  routine  for  the  inverse 
problem  to  determine  the  endpoints  of  a  uniform  distribution  of  relaxation  times  for  each  set  of  initial 
conditions  (case  0  represents  the  exact  solution).  The  values  in  parenthesis  denote  the  absolute  value  of  the 
difference  as  the  number  of  digits  shown  here  would  not  suffuciently  distinguish  the  approximations  from 
the  exact  solution.  _ 


case 

1011 

Frequency  (Hz) 
109 

106 

0 

0.248369 

0.248369 

5.42467e-8 

1 

(1.64144e-8) 

0.248701 

5.43479e-8 

2 

(1.51187e-9) 

0.248396 

5.43315e-8 

3 

(5.12895e-ll) 

0.248373 

5.42813e-8 

Table  48:  Resulting  values  of  ra  from  minimizing  along  the  line  of  constant  A  (or  r  for  /  =  106  Hz),  for  the 
inverse  problem  to  determine  the  endpoints  of  a  uniform  distribution  of  relaxation  times  (recall  the  exact 
solution  r*  =3.16000e-8). 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

3.08694e-08 

3.1653e-08 

3.17617e-08 

2 

3.16401e-08 

3.16043e-08 

3.17357e-08 

3 

3.15905e-08 

3.16007e-08 

3.16559e-08 

Finally,  for  fine  tuning,  we  apply  the  full  two  parameter  Levenberg-Marquardt  routine  using  the  estimates 
from  minimizing  along  constant  A  (or  f).  These  results  are  shown  in  Tables  50  and  51.  We  see  that  the 
estimates  change  very  little,  if  at  all,  which  suggests  that  our  approximation  method  is  not  only  efficient, 
but  quite  accurate  as  well. 
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Table  49:  Resulting  values  of  t>,  from  the  minimizing  along  the  line  of  constant  A  (or  f  for  /  =  10 6Hz), 
for  the  inverse  problem  to  determine  the  endpoints  of  a  uniform  distribution  of  relaxation  times  (recall  the 
exact  solution  rb  =1.5800e-8). 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

1.68829e-08 

1.56448e-08 

1.55282e-08 

2 

1.57433e-08 

1.57874e-08 

1.55715e-08 

3 

1.58134e-08 

1.5798e-08 

1.57053e-08 

Table  50:  Resulting  values  of  r0  from  the  two  parameter  Levenberg-Marquardt  routine  for  the  inverse 
problem  to  determine  the  endpoints  of  a  uniform  distribution  of  relaxation  times  (recall  the  exact  solution 
r*  =3.16000e-8).  _ 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

3.08694e-08 

3.15996e-08 

3.16014e-08 

2 

3.16401e-08 

3.16005e-08 

3.16008e-08 

3 

3.15905e-08 

3.16008e-08 

3.16002e-08 

Table  51:  Resulting  values  of  rb  from  the  two  parameter  Levenberg-Marquardt  routine  for  the  inverse 
problem  to  determine  the  endpoints  of  a  uniform  distribution  of  relaxation  times  (recall  the  exact  solution 
rb*  =1.58000e-8).  _ 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

1.68829e-08 

1.58005e-08 

1.57988e-08 

2 

1.57433e-08 

1.57993e-08 

1.57991e-08 

3 

1.58134e-08 

1.57988e-08 

1.57998e-08 

7  Inverse  Problem  Results  Using  a  Gaussian  Distribution 


In  this  section  we  explore  the  possibility  of  determining  dielectric  parameters  that  are  normally  distributed. 
We  apply  a  Gaussian  distribution  to  the  logarithms  of  the  relaxation  times  (thus  actually  the  relaxation  times 
are  log-normally  distributed,  but  all  of  our  computions  are  in  “log”-space  so  we  refer  to  these  parameters  as 
normally  distributed).  Based  on  experimental  data  in  [23],  this  seems  to  be  a  likely  model  for  the  relaxation 
times  of  a  material. 

Recall  from  Section  3.3  that  our  PDF  is 


dF(r;  /x,  a) 


1  11 
\/2tt(t2  ln  10  t 


exp 


(log 

— 2, ? - )dT 


(24) 


for  a  log-normal  distribution  in  r.  Here  /z  is  the  mean  of  the  relaxation  times  and  a  is  the  standard  deviation. 
Recall  further  that  the  density  is  truncated  to  a  support  interval  [ra,rb]  where  ra  and  rb  are  determined 
based  on  /z  and  er  as  follows 


Ta  :=  10M_6<T 
n  :=  10At+6CT. 

Thus,  this  problem  again  presents  us  with  a  two  parameter  inverse  problem,  namely,  to  try  to  resolve  the 
mean  and  standard  deviation  of  the  distribution  used  to  generate  the  given  data.  Learning  from  the  success 
in  the  previous  section,  we  expect  to  apply  a  three  step  approach:  one  parameter  Levenberg-Marquardt  for 
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the  mean  p  (which  we  expect  to  converge  to  the  line  of  constant  A  or  f),  then  traverse  the  line  of  constant 
A  or  f  (which  is  again  one  parameter),  and  finally,  perform  two  parameter  Levenberg-Marquardt  for  fine- 
tuning.  We  can  see  from  the  surface  plot  using  /  =  1011  shown  in  Figure  18  that  we  do  again  have  a  trench 
representing  a  “line  of  best  fit”,  which  does  in  fact  have  a  single  minimizer,  but  it  appears  to  be  independent 
of  a.  The  trench  is  actually  slightly  slanted,  as  can  be  seen  in  Figure  19  which  zooms  in  on  /i,  and  in  fact,  for 
the  lower  frequencies  the  slant  is  indeed  in  the  opposite  direction.  Therefore  while  at  first  glance  of  Figure 
18  traversing  a  line  of  constant  A  or  t  may  seem  unnecessary,  it  does  actually  give  better  results  than  simply 
minimizing  in  the  a  direction  alone. 

Remark  2  Note  that  the  slant  does  not  appear  to  be  a  consequence  of  our  quadrature  rule.  Changing  the 
number  of  nodes  used  in  our  computations  had  no  effect  on  the  slope.  Further,  increasing  the  support  on 
which  the  integral  was  computed  had  no  effect,  e.g.,  to  12a  on  either  side  of  p.  While  f  Q(0,z)  =  p  is 
actually  independent  of  a  we  still  expect  a  to  affect  the  objective  function  in  some  other  way,  and  this  is  in 
fact  what  is  happening.  Lastly,  the  apparent  local  minima  in  Figure  19  are  again  just  an  artifact  of  the  grid 
points  not  falling  exactly  on  the  “line  of  best  fit”. 

We  performed  our  inverse  problem  solution  approach  to  the  sample  problem  of  p  =  xo-' 62525  and 
a  =  0.0457575  (this  choice  of  parameters  results  in  a  distribution  function  comparable  to  the  uniform 
distribution  case  of  the  previous  section  in  that  the  corresponding  densities  have  the  same  mean  and  roughly 
similar  support,  see  Figure  20).  The  cases  of  initial  values  we  considered  for  p  and  a,  and  the  corresponding 
initial  A  and  f  values,  are  given  in  Table  52.  Note  that  Case  3  is  specifically  designed  to  test  whether  the 
standard  deviation  can  be  determined  when  the  mean  is  known. 


Table  52:  Initial  estimates  and  corresponding  A  and  t  values  for  the  inverse  problem  to  determine  the  mean 
and  standard  deviation  of  a  normal  distribution  of  relaxation  times  (case  0  corresponds  to  the  exact  solution) . 


case 

log  (/r0) 

co 

Ao 

0 

-7.62525 

0.0457575 

0.141524 

2.38319e-08 

1 

-6.92628 

0.036606 

0.0282483 

1.18922e-07 

2 

-8.32422 

0.054909 

0.709351 

4.77804e-09 

3 

-7.62525 

0.0411817 

0.141375 

2.37329e-08 

The  p  estimates  resulting  from  the  one  parameter  Levenberg-Marquardt  are  given  in  Table  53.  In  each 
case  at  least  three  decimals  places  of  accuracy  was  achieved.  Also,  at  least  seven  decimal  places  of  accuracy 
are  obtained  on  the  estimates  of  the  corresponding  A  and  f  values,  as  displayed  in  Table  54.  The  p  and  a 
estimates  resulting  from  the  one  parameter  tracing  of  the  “line  of  best  fit”  are  given  in  Tables  55  and  56, 
respectively.  In  each  case  at  least  six  decimals  places  of  accuracy  was  achieved  for  the  log  of  p.  Also,  at  least 
three  decimal  places  of  accuracy  are  obtained  on  the  estimates  of  a.  Finally,  the  p  and  er  estimates  resulting 
from  the  two  parameter  Levenberg-Marquardt  fine-tuning  are  given  in  Tables  57  and  58,  respectively.  In  all 
but  a  couple  cases  we  have  as  much  degree  accuracy  as  the  number  of  digits  shown  here  allows.  Figure  21 
gives  a  graphical  representation  of  how  well  our  inverse  problem  solution  method  converged  given  the  initial 
condition  Case  1  from  Table  52.  We  conclude  that  the  inverse  problem  involving  a  Gaussian  distribution 
of  relaxation  times  for  a  Debye  polarization  model  is  computationally  feasible  for  the  sample  parameters 
presented  here. 
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Log  of  objective  function  for  gaussian  distribution  using  f=1e+1 1 
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Figure  18:  The  log  of  the  objective  function  for  the  gaussian  distribution  inverse  problem  versus  the  log  of 
and  a  using  a  frequency  of  lO^Hz. 

log  of  objective  function  for  f=1e+1 1 


Figure  19:  The  log  of  the  objective  function  for  the  gaussian  distribution  inverse  problem  versus  the  log  of 
/i  and  a,  zoomed  in  with  respect  to  /i,  and  using  a  frequency  of  lO^Hz.  The  appearance  of  several  local 
mininima  is  merely  an  artifact  of  the  lattice  points  used  to  plot  the  surface. 
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U(xa=1 .58e-08,xb=3.16e-08)  and  LN(p  =-7.62525, a  =0.0457575) 
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Figure  20:  This  plot  compares  the  uniform  density  considered  in  Section  6  to  the  lognormal  density.  Shown 
are  the  two  density  functions  and  the  corresponding  quadrature  nodes  used  in  integrating  them. 


Table  53:  The  resulting  log(//)  values  from  performing  a  one  parameter  Levenberg-Marquardt  optimization 
for  the  inverse  problem  to  determine  the  mean  and  standard  deviation  of  a  normal  distribution  of  relaxation 
times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is  log(/i*)  =  —7.62525). 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

-7.62612 

-7.6261 

-7.6244 

2 

-7.62419 

-7.62422 

-7.62629 

3 

-7.62571 

-7.62701 

-7.62348 

Table  54:  The  corresponding  A  and  f  values  from  performing  a  one  parameter  Levenberg-Marquardt  op¬ 
timization  for  the  inverse  problem  to  determine  the  mean  and  standard  deviation  of  a  normal  distribution 
of  relaxation  times  for  each  frequency  and  initial  value  case.  Note  that  the  absolute  values  of  the  differ¬ 
ences  between  the  estimates  and  the  exact  values  are  shown  as  the  number  of  digits  shown  here  would  not 
sufficiently  distinguish  the  approximations  from  the  exact  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

2.62469e-ll 

7.38776e-08 

8.98789e-13 

2 

3.61986e-ll 

9.10257e-08 

1.10726e-12 

3 

4.25993e-13 

1.56168e-07 

1.86751e-12 

38 


Table  55:  The  resulting  log(^)  values  from  traversing  the  “line  of  best  fit”  for  the  inverse  problem  to 
determine  the  mean  and  standard  deviation  of  a  normal  distribution  of  relaxation  times  for  each  frequency 
and  initial  value  case  (recall  the  exact  solution  is  log(/i*)  =  —7.62525). 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

-7.62526 

-7.62525 

-7.62526 

2 

-7.62524 

-7.62525 

-7.62524 

3 

-7.62524 

-7.62525 

-7.62527 

Table  56:  The  resulting  a  values  from  traversing  the  “line  of  best  fit”  for  the  inverse  problem  to  determine 
the  mean  and  standard  deviation  of  a  normal  distribution  of  relaxation  times  for  each  frequency  and  initial 
value  case  (recall  the  exact  solution  is  a*  =  0.0457575). 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

0.0456376 

0.0457564 

0.0456871 

2 

0.0458453 

0.0457592 

0.0458442 

3 

0.0458207 

0.0457551 

0.0456102 

Table  57:  The  resulting  log  (fi)  values  from  performing  a  two  parameter  Levenberg-Marquardt  optimization 
for  fine-tuning  of  the  solution  for  the  inverse  problem  to  determine  the  mean  and  standard  deviation  of  a 
normal  distribution  of  relaxation  times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is 
log  (/u*)  =  -7.62525). 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

-7.62526 

-7.62525 

-7.62525 

2 

-7.62524 

-7.62525 

-7.62525 

3 

-7.62524 

-7.62525 

-7.62525 

Table  58:  The  resulting  a  values  from  performing  a  two  parameter  Levenberg-Marquardt  optimization  for 
fine-tuning  of  the  solution  for  the  inverse  problem  to  determine  the  mean  and  standard  deviation  of  a 
normal  distribution  of  relaxation  times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is 
cr*  =  0.0457575). 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

0.0456376 

0.0457576 

0.0457575 

2 

0.0458453 

0.0457576 

0.0457575 

3 

0.0458207 

0.0457575 

0.0457575 
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Estimated  density  of  x  as  log  normal 


Figure  21:  Shown  are  the  initial  density  function,  the  minimizing  density  function  and  the  true  density 
function  (the  latter  two  being  practically  identical)  corresponding  to  Case  1  in  Table  52. 


40 


8  Inverse  Problem  Results  Using  a  Bi-Gaussian  Distribution 

In  this  section  we  combine  the  two  ideas  from  Sections  5  and  7  in  that  we  consider  a  material  comprised  of 
primarily  two  distinct  polarization  mechanisms,  but  rather  than  assuming  two  atoms  at  t\  and  r2  as  before 
we  let  /ii  =  7”i  and  /i2  =  t2  in  two  Gaussian  distributions.  Essentially,  we  are  decomposing  the  distribution 
function  into  two  components,  each  dependent  on  distinct  means  and  standard  deviations  as  follows: 

dF  =  adF(r ;  fj,1}  cy)  +  (1  -  a)dF(r;  /z2,  er2),  (25) 

where  dF  is  given  by  (24)  (we  take  a  =  .5  in  our  numerical  experiments). 

We  now  have  a  four  parameter  inverse  problem,  namely 

min  J(F) 

where  Q  is  the  admissible  region  for  our  unknown  parameters  (e.g.,  a i  >0).  For  the  following  examples,  the 
true  values  of  the  means  were  taken  to  be  (/r*,/r2)  =  (10— 7'80134, 10  ' ■50031),  with  the  standard  deviations 
ranging  from  .02  to  .07  As  in  the  previous  gaussian  case,  we  expect  the  objective  function  to  be  relatively 
less  dependent  on  the  standard  deviations  than  the  means,  therefore  we  address  the  dependency  on  the 
means  first.  Figure  22  shows  a  surface  plot  of  the  objective  function  with  respect  to  /ii  and  /i2.  Here  we 
can  see  the  similarities  to  the  distribution  function  from  Section  5,  namely  the  presense  of  the  “line  of  best 
fit”  with  two  (symmetric)  global  minima.  (When  comparing  to  Figure  3,  note  that  here  we  have  rotated  the 
axis  for  better  viewing  of  the  minima  under  the  surface.)  Our  minimization  approach  is  thus  the  same  as 
that  which  worked  well  in  the  previous  sections.  We  perform  a  one  parameter  Levenberg-Marquardt  search 
in  the  /ii  direction  (and  just  to  be  sure,  a  one  parameter  Levenberg-Marquardt  search  in  the  /r2  direction), 
then  optimize  along  the  “line  of  best  fit” ,  and  finally  fine-tune  with  a  two  parameter  Levenberg-Marquardt 
search  for  both  /ri  and  /i2. 


log  of  objective  function  for  bigaussian  distribution  using  f=1  e+1 1 


Figure  22:  The  log  of  the  objective  function  for  the  bi-gaussian  distribution  inverse  problem  versus  the  log 
of  /j,i  and  H2  using  a  frequency  of  10 11  Hz. 


41 


Thus  far  we  have  completely  ignored  the  standard  deviations.  Figures  23  and  24  show  the  objective 
function  versus  /ri  and  /i2  for  two  distinct  cases:  using  the  exact  values  for  the  standard  deviations  and 
using  significantly  incorrect  ones.  The  fact  that  the  location  of  the  “line  of  best  fit”  does  not  change 
drastically  suggests  that  we  should  be  able  to  at  least  determine  the  values  of  A  or  f  with  just  our  one 
parameter  searches,  as  in  the  previous  sections.  However,  the  location  of  the  local  minima  along  this  curve 
has  changed,  therefore  we  should  not  expect  our  /xi  and  /r 2  estimates  so  far  to  be  our  final  answer.  We 
must  first  attempt  to  optimize  with  respect  to  the  standard  deviations,  and  then  finally  attempt  a  full  four 
parameter  minimization  to  make  sure  we  have  accounted  for  all  interdependencies.  These  are  the  last  two 
steps  of  our  now  six  step  optimization  process  for  the  four  parameter  inverse  problem  involving  a  Bi-Gaussian 
distribution. 


Figure  23:  The  log  of  the  objective  function  for  the  bi-gaussian  distribution  inverse  problem  versus  the  log 
of  /ri  and  /Z2  using  correct  values  for  the  standard  deviations. 


log  of  objective  function  for  bigaussian  distribution  using  f=1e+1 1  with  incorrect  a 
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Figure  24:  The  log  of  the  objective  function  for  the  bi-gaussian  distribution  inverse  problem  versus  the  log 
of  /j,i  and  112  using  incorrect  values  for  the  standard  deviations. 


42 


The  true  values  of  the  optimization  parameters  and  the  initial  conditions  considered  for  each  case  we 
attempted  are  given  in  Table  59.  These  values  correspond  to  roughly  the  same  relative  initial  error  as  the 
initial  value  cases  in  the  single  gaussian  problem  of  Section  7.  The  pi  and  /i2  estimates  resulting  from  the 
one  parameter  Levenberg-Marquardt  searches  are  given  in  Tables  60  and  61,  respectively. 

Remark  3  Note  that  in  Case  2  at  f  =  1011  the  value  of  p±  diverged  toward  zero  in  our  first  attempts. 
Therefore  the  residts  shown  in  the  tables  for  Case  2  at  each  frequency  represent  first  minimizing  in  the  /z2 
direction  followed  by  minimization  in  the  p±  direction.  In  general  better  results  are  achieved  when  optimizing 
the  smaller  of  the  values  (i.e.,  corresponding  to  a  more  negative  logarithm)  first. 

Several  of  the  /i2  estimates  did  not  change  from  their  initial  values  during  the  one  parameter  optimization 
for  /i2  suggesting  that  the  p\  search  arrived  sufficiently  close  to  the  A  or  f  curve  (likewise  for  Case  2,  the 
two  higher  frequency  estimates  for  pi  did  not  change  suggesting  that  the  /i2  search  converged  sufficiently). 
The  values  of  A  or  t  from  the  two  one-parameter  searches  are  given  in  Table  62.  These  values,  on  average, 
have  about  two  fewer  decimal  places  of  accuracy  as  those  in  the  single  gaussian  case.  Still,  as  the  following 
results  shall  show,  in  general  this  is  sufficiently  close  to  allow  for  good  final  estimates. 

We  now  hold  A  or  t  fixed  and  search  for  p\  and  /i2.  The  p \  and  /i2  estimates  resulting  from  the  one 
parameter  tracing  of  the  “line  of  best  fit”  are  given  in  Tables  63  and  64,  respectively.  Finally,  the  p±  and  p2 
estimates  resulting  from  the  two  parameter  Levenberg-Marquardt  fine-tuning  are  given  in  Tables  65  and  66, 
respectively.  The  lack  of  a  significant  difference  after  the  two  parameter  search  suggests  that  we  have  found 
local  minima  with  respect  to  pi  and  /i2  given  our  initial  standard  deviation  estimates. 

Remark  4  Note  that  in  some  cases  (actually,  about  half!)  p\  is  converging  toward  (we  have  highlighted 
these  cases  in  each  subsequenct  table).  This  is  the  same  problem  that  we  encountered  in  the  discrete  case 
in  Section  5.  For  our  examples  here,  a  =  .5  so  the  problem  is  symmetric  provided  the  standard  deviations 
converge  to  the  corresponding  symmetric  solution  as  well  (i.e.,  a®  converges  to  a \).  We  could  restrict  our 
optimization  problem  to,  for  instance,  only  allow  p±  >  P2,  but  we  perfer  the  idea  that  there  are  two  equally 
viable  solutions  to  choose  from  since  it  may  actually  double  our  chances  of  finding  one  of  them!  If  our  initial 
estimates  for  each  a were  based  on  a  certain  ordering  of  the  pt ’s  (for  example,  say  we  know  the  distribution 
of  the  smaller  relaxation  time  is  much  more  narrow)  we  could  easily  test  for  that  ordering  and  rearrange  our 
g i  estimates  if  necessary  (though  we  would  need  to  optimize  with  respect  to  each  pi  once  more  to  reflect  the 
changes). 

While  the  results  for  determining  the  means  of  the  unknown  Bi-Gaussian  distribution  are  surprisingly 
decent,  the  optimization  of  the  standard  deviations  was  expectedly  more  difficult.  Figure  25  demonstrates 
that  the  objective  function  behaves  in  a  similar  way  with  respect  to  the  standard  deviations  as  it  did  to  the 
means  in  that  there  is  a  “line  of  best  fit”.  Unfortunately,  when  compared  to  Figure  22  it  becomes  obvious 
that  the  scale  involved  in  traversing  along  this  curve  is  too  small  to  be  accurate  enough  for  optimization. 
This  is  in  fact  what  we  experienced  in  attempting  to  minimize  with  respect  to  g\  along  the  “line  of  best 
fit” .  The  changes  in  each  a±  were  in  almost  all  cases  less  than  six  decimal  places.  Applying  a  two  parameter 
Levenberg-Marquardt  for  fine-tuning  afterwards  gave  slightly  better  results  for  some  cases,  and  we  display 
these  approximations  in  Tables  67  and  68.  Using  these  values  along  with  the  best  estimates  so  far  for  our 
means,  we  run  a  four  parameter  Levenberg-Marquardt  search.  The  results  for  this  procedure  are  given  in 
Tables  69  through  72. 

One  may  argue  that  the  lack  of  sensitivity  to  a  may  be  due  to  starting  initial  estimates  too  close  to 
the  true  solution,  but  this  would  not  explain  why,  with  all  other  parameters  held  fixed,  gradient  based 
methods  take  the  standard  deviation  estimates  farther  from  the  true  solution  in  half  of  the  cases  we  tried. 
It  is  more  likely  that  the  numerical  errors  involved  in  the  simulations  are  affecting  the  gradients.  However, 
the  results  here  reflect  attempts  that  have  been  made  to  account  for  this  using  varying  gradient  step  sizes 
(our  Levenberg-Marquardt  routine  exits  only  after  verifying  small  gradients  with  three  different  step  sizes 
as  explained  in  [13]). 

Our  conclusion  for  the  problem  of  determining  the  means  and  the  standard  deviations  for  a  bi-gaussian 
distribution  of  relaxation  times  is  that  the  means  are  very  readily  determined  with  reasonable  initial  esti¬ 
mates,  even  with  significantly  incorrect  standard  deviations.  However,  the  standard  deviations  are  not  easily 
determined  even  with  quite  accurate  estimates  for  the  means. 
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Table  59:  Initial  estimates  for  the  inverse  problem  to  determine  the  means  and  standard  deviations  of  a 
bi-gaussian  distribution  of  relaxation  times.  For  each  case,  (/z*,  /z£)  =  (10-7  80134,  lO-7  '50031).  For  Case  1, 
=  (0.0457575,0.0457575).  For  Cases  2  and  3,  (oj^oj)  =  (0.0705811,0.0222764). 


case 

log(Mi) 

log(M2) 

02 

i 

-8.50031 

0.036606 

-6.80134 

0.0571969 

2 

-7.10237 

0.0846973 

-8.19928 

0.0185637 

3 

-8.50031 

0.0352905 

-6.80134 

0.0445528 

Table  60:  The  resulting  log(/zi)  values  from  performing  a  one  parameter  Levenberg-Marquardt  optimization 
in  Hi  for  the  inverse  problem  to  determine  the  means  and  standard  deviations  of  a  bi-gaussian  distribution 
of  relaxation  times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is  log(/zJ)  =  —7.80134). 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

-7.94848 

-7.94847 

-8.37177 

2 

-7.10237 

-7.10237 

-7.39253 

3 

-7.95035 

-7.95034 

-8.37264 

Table  61:  The  resulting  log(/i2)  values  from  performing  a  one  parameter  Levenberg-Marquardt  optimization 
in  H2  for  the  inverse  problem  to  determine  the  means  and  standard  deviations  of  a  bi-gaussian  distribution 
of  relaxation  times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is  log(/iJ)  =  —7.50031). 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

-6.80134 

-6.80377 

-7.36944 

2 

-7.91764 

-7.91763 

-8.13781 

3 

-6.80134 

-6.80375 

-7.36803 

Table  62:  The  corresponding  A  and  t  values  from  performing  each  one  parameter  Levenberg-Marquardt 
optimization  in  /zi  and  M2  for  the  inverse  problem  to  determine  the  means  and  standard  deviations  of  a 
bi-gaussian  distribution  of  relaxation  times  for  each  frequency  and  initial  value  case.  Note  that  the  absolute 
values  of  the  differences  between  the  estimates  and  the  exact  values  arc  shown. 


case 

1011 

Frequency  (Hz) 
109 

106 

i 

7.945e-10 

5.7109e-05 

1.45861e-10 

2 

5.00444e-10 

1.73813e-06 

4.24458e-10 

3 

1.74996e-08 

5.57019e-05 

1.46209e-10 
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Table  63:  The  resulting  log(/ri)  values  from  traversing  the  “line  of  best  fit”  for  the  inverse  problem  to 
determine  the  mean  and  standard  deviation  of  a  normal  distribution  of  relaxation  times  for  each  frequency 
and  initial  value  case  (recall  the  exact  solution  is  log(/U*)  =  —7.80134).  The  highlighted  values  denote 
convergence  to  the  symmetric  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

1 

-7.80682 

-7.5052 

-7.79549 

2 

-7.47853 

-7.47881 

-7.50275 

3 

-7.81476 

-7.48674 

-7.79517 

Table  64:  The  resulting  log(/i2)  values  from  traversing  the  “line  of  best  fit”  for  the  inverse  problem  to 
determine  the  mean  and  standard  deviation  of  a  normal  distribution  of  relaxation  times  for  each  frequency 
and  initial  value  case  (recall  the  exact  solution  is  log^)  =  —7.50031).  The  highlighted  values  denote 
convergence  to  the  symmetric  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

1 

-7.48957 

-7.79775 

-7.50821 

2 

-7.81299 

-7.81285 

-7.78246 

3 

-7.47935 

-7.8106 

-7.50706 

Table  65:  The  resulting  log(/ii)  values  from  performing  a  two  parameter  Levenberg-Marquardt  optimization 
for  fine-tuning  of  (^1,^2)  for  the  inverse  problem  to  determine  the  means  and  standard  deviations  of  a  bi- 
gaussian  distribution  of  relaxation  times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution 
is  log(/z*)  =  —7.80134).  The  highlighted  values  denote  convergence  to  the  symmetric  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

1 

-7.80682 

-7.5052 

-7.79549 

2 

-7.47851 

-7.47881 

-7.50275 

3 

-7.81442 

-7.48674 

-7.79517 

Table  66:  The  resulting  log(/i2)  values  from  performing  a  two  parameter  Levenberg-Marquardt  optimization 
for  fine-tuning  of  ( fii,l-i2 )  for  the  inverse  problem  to  determine  the  means  and  standard  deviations  of  a  bi- 
gaussian  distribution  of  relaxation  times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution 
is  log^)  =  —7.50031).  The  highlighted  values  denote  convergence  to  the  symmetric  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

1 

-7.48957 

-7.79775 

-7.50821 

2 

-7.813 

-7.81285 

-7.78246 

3 

-7.47962 

-7.8106 

-7.50706 
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Figure  25:  The  log  of  the  objective  function  for  the  bi-gaussian  distribution  inverse  problem  versus  <J\  and 
a 2  using  a  frequency  of  1011  Hz. 


Table  67:  The  resulting  o\  values  from  a  two  parameter  Levenberg-Marquardt  search  for  g\  and  a-2  for  the 
inverse  problem  to  determine  the  means  and  standard  deviations  of  a  bi-gaussian  distribution  of  relaxation 
times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is  g\  =  0.0457575  for  Case  1  and 
0.0705811  otherwise).  The  highlighted  values  denote  convergence  to  the  symmetric  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

1 

0.036606 

0.00603311 

0.0136791 

2 

0.0846974 

0.0844293 

0.0761876 

3 

0.0352905 

0.00857058 

0.0374784 

Table  68:  The  resulting  02  values  from  a  two  parameter  Levenberg-Marquardt  search  for  g\  and  02  for  the 
inverse  problem  to  determine  the  means  and  standard  deviations  of  a  bi-gaussian  distribution  of  relaxation 
times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is  g%  =  0.0457575  for  Case  1  and 
0.0222764  otherwise).  The  highlighted  values  denote  convergence  to  the  symmetric  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

1 

0.0571969 

0.060935 

0.0658016 

2 

0.0185637 

0.0192023 

0.044288 

3 

0.0445528 

0.0482781 

0.0456734 

46 


Table  69:  The  resulting  log(/zi)  values  from  performing  a  four  parameter  Levenberg-Marquardt  optimization 
for  all  parameters  in  the  inverse  problem  to  determine  the  means  and  standard  deviations  of  a  bi-gaussian 
distribution  of  relaxation  times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is  log(/x^)  = 
—7.80134).  The  highlighted  values  denote  convergence  to  the  symmetric  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

1 

-7.80682 

-7.5052 

-7.78543 

2 

-7.47851 

-7.47911 

-7.50072 

3 

-7.81432 

-7.48831 

-7.78507 

Table  70:  The  resulting  \og(/i2)  values  from  performing  a  four  parameter  Levenberg-Marquardt  optimization 
for  all  parameters  in  the  inverse  problem  to  determine  the  means  and  standard  deviations  of  a  bi-gaussian 
distribution  of  relaxation  times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is  log^O)  = 
—7.50031).  The  highlighted  values  denote  convergence  to  the  symmetric  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

1 

-7.48957 

-7.79687 

-7.50993 

2 

-7.813 

-7.81279 

-7.80691 

3 

-7.47506 

-7.80992 

-7.50848 

Table  71:  The  resulting  u\  values  from  a  four  parameter  Levenberg-Marquardt  search  for  all  parameters 
in  the  inverse  problem  to  determine  the  means  and  standard  deviations  of  a  bi-gaussian  distribution  of 
relaxation  times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is  a  l  =  0.0457575  for 
Case  1  and  0.0705811  otherwise).  The  highlighted  values  denote  convergence  to  the  symmetric  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

1 

0.0366059 

0.0056312 

0.0136811 

2 

0.0846973 

0.0832475 

0.0736677 

3 

0.0352905 

0.00879725 

0.037513 

Table  72:  The  resulting  02  values  from  a  four  parameter  Levenberg-Marquardt  search  for  all  parameters 
in  the  inverse  problem  to  determine  the  means  and  standard  deviations  of  a  bi-gaussian  distribution  of 
relaxation  times  for  each  frequency  and  initial  value  case  (recall  the  exact  solution  is  cr|  =  0.0457575  for 
Case  1  and  0.0222764  otherwise).  The  highlighted  values  denote  convergence  to  the  symmetric  solution. 


case 

1011 

Frequency  (Hz) 
109 

106 

1 

0.0571968 

0.0604793 

0.0663628 

2 

0.0185637 

0.0198981 

0.0442772 

3 

0.0445528 

0.0477348 

0.0463262 
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9  Summary 

In  this  report  we  have  discussed  theoretical  and  numerical  results  for  inverse  problems  involving  Maxwell’s 
equations  with  a  general  polarization  term  which  includes  probability  distributions  for  both  dielectric  pa¬ 
rameters  and  mechanisms.  A  theoretical  framework  to  treat  such  classes  of  problems  was  given  along  with 
stability  results  for  the  corresponding  least  squares  inverse  problems.  We  presented  examples  of  distribu¬ 
tions  of  parameters  including  discrete,  uniform,  log-normal,  and  log-bi-gaussian.  In  each  case  the  objective 
function  was  characterized  by  a  “line  of  best  fit”  corresponding  to  a  curve  of  constant  value  of  either  A  or  r 
depending  on  the  relative  values  of  the  interrogating  frequency  ui,  and  1/f,  where  f  is  the  weighted  average 
of  the  relaxation  times.  The  constant  values  of  these  curves  turned  out  to  be,  in  fact,  t  and  A,  which  is  the 
weighted  average  of  the  inverses  of  the  relaxation  times  (scaled  by  the  speed  of  light  for  convenience) .  The 
similarities  between  the  cases  were  exploited  in  determining  effective  optimization  procedures  for  subsequent 
cases. 

For  the  discrete  case  involving  two  atoms,  each  value  of  r  was  readily  determined  if  the  volume  propor¬ 
tions  were  known.  Likewise,  the  volume  proportion  could  be  determined  given  the  values  of  each  relaxation 
time.  However,  the  problem  of  determining  all  three  parameters  simultaneously  is  under-determined.  Al¬ 
though,  in  certain  situations,  interrogating  with  a  different  frequency  (in  particular,  greater  than  the  critical 
frequency  if  the  original  interrogating  frequency  is  smaller)  can  provide  enough  new  information  for  the 
volume  proportions  and  the  relaxation  times  to  be  determined  simultaneously. 

The  uniform  distribution  is  generally  straight-forward  to  optimize;  each  endpoint  of  the  distribution  was 
determined,  on  average,  to  about  three  decimal  places  of  accuracy. 

The  inverse  problem  involving  the  log-normal  distribution,  or  rather  a  gaussian  distribution  on  the 
logarithms  of  the  relaxation  times,  behaved  similarly  to  the  discrete  distribution,  and  therefore  the  previous 
solution  methods  worked  well.  The  mean  of  the  distribution  was  determined  to  about  four  decimal  places  in 
log-space.  The  standard  deviation  was  also  determined  quite  accurately  (about  three  decimal  places),  but 
it  was  evident  in  the  surface  plots  that  the  objective  function  was  particularly  insensitive  to  the  standard 
deviation.  Any  amount  of  noise  in  the  system  may  prevent  accurate  results. 

Finally,  we  considered  a  log-bi-gaussian  distribution  of  relaxation  times.  As  in  the  gaussian  distribution 
problem,  the  objective  function  was  insensitive  to  the  standard  deviations.  Therefore,  even  with  quite 
accurate  estimates  for  the  means,  no  usable  estimates  for  the  standard  deviations  could  be  found.  However, 
in  contrast,  rather  accurate  solutions  for  both  of  the  means  could  be  determined  using  substantially  incorrect 
estimates  for  the  standard  deviations.  Still,  the  effect  of  the  standard  deviation  on  the  objective  function 
is  sufficient  to  be  able  to  distinguish  distributions  which  are  fairly  broad  (large  standard  deviations)  from 
those  that  are  nearly  delta  functions  (very  small  standard  deviations).  This  is  important  as  it  suggests  that 
systems  which  truly  have  continuous  distributions  for  their  parameters  should  not  be  modeled  using  discrete 
distributions.  Further,  an  objective  function  of  the  type  in  this  report  may  be  used  to  create  a  measure  of 
the  error  induced  in  the  signal  by  opting  to  use  a  discrete  distribution  instead  of  a  continuous  one  when  it 
is  sufficiently  helpful  for  computational  simplicity  to  do  so. 

In  fact,  homogenization  techniques  determine  single  values  of  relaxation  times  which  produce  signals  that 
closely  match  those  generated  with  more  complex  mechanisms  (such  as  a  distribution  of  Debye  polarization 
mechanisms)-for  example,  see  [7].  It  turns  out  that  from  our  analysis  here  that  for  lower  frequencies  {u>l  < 
ujc  =  1/f )  the  effective  relaxation  time  is  simply  the  weighted  average  of  the  actual  relaxation  times  (or  for 
continuous  distributions,  the  integral  of  r  with  respect  to  the  distribution).  However,  for  higher  frequencies 
(u>h  >  ojc  —  1/t)  the  effective  relaxation  time  is  actually  the  inverse  of  the  weighted  average  of  inverses  (or 
the  inverse  of  the  integral  of  1/r  with  respect  to  the  distribution  in  the  continuous  case).  For  example, 

1 


where  a*  represents  the  volume  fraction  of  the  material  with  discrete  relaxation  times  7y.  Both  the  low 
frequency  and  high  frequency  estimates  are  common  results  in  homogenization  theory.  However,  any  ho¬ 
mogenization  procedure  which  does  not  take  the  interrogating  frequency  into  account  will  not  be  able  to 
determine  which  is  the  correct  effective  relaxation  time  for  the  situation. 

Lastly,  it  should  be  noted  that  the  procedure  of  traversing  a  “line  of  best  fit”,  as  done  in  this  report,  is 
only  a  valid  approach  to  solving  non-linear  optimization  problems  (which  have  relations  that  are  relatively 
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optimal)  if  the  equation  of  this  curve  is  known.  Other  possible  solutions  methods  for  more  general  problems 
include  taking  an  “arbitrary  orthogonal  step” ,  or  resorting  to  a  simplex  search  for  those  parts  of  the  inverse 
problem  that  are  particularly  difficult  for  the  gradient  based  methods. 
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